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Abstract.  Energy  compaction  filters  have  attracted  considerable  attention  due  in  part,  to  the  fact  that  they 
are  the  building  blocks  of  optimal  orthonormal  (paraunitary)  filter  banks.  In  this  paper  we  introduce  some  new 
design  techniques  for  optimum  M-channel  FIR  compaction  filters  for  a  given  input  power  spectrum.  Some 
properties  of  the  optimum  FIR  compaction  filters  and  the  corresponding  maximum  compaction  gains  are  also 
derived.  For  the  design  part,  a  modification  of  the  well-known  linear  programming  technique  is  considered. 
We  also  consider  multistage  (IFIR)  designs  of  compaction  filters.  A  new,  efficient  design  method  called  the 
window  method  is  then  introduced.  The  method  generates  M— channel  FIR  compaction  filters  for  an}'  given 
power  spectrum.  Although  it  is  suboptimal,  no  optimization  tools  of  any  kind  are  involved  and  the  algorithm 
terminates  in  a  finite  number  of  elementary  steps.  As  the  filter  order  increases,  the  window  method  produces 
compaction  gains  that  are  very  close  to  the  optimal  ones.  We  give  a  necessary  condition  for  a  compaction  filter 
to  be  optimum  and  provide  some  bounds  on  the  maximum  compaction  gains.  Finally  we  propose  an  analyical 
method  for  the  two-channel  case  which  finds  the  optimum  FIR  compaction  filters  for  a  class  of  random  processes. 


DTIC  QUALITY  INSPECTED  2 


t  Work  supported  in  parts  by  Office  of  Naval  Research  grant  N00014-93-1-0231,  Tektronix,  Inc.,  and  Rockwell 
Inti. 

DISTRIBUTION  STATEMENT"~A 

Approved  for  public  release; 

Distribution  Unlimited 


I.  INTRODUCTION 

Energy  compaction  filters  have  attracted  a  great  deal  of  attention  due  in  part,  to  the  fact  that  they  are 
the  building  blocks  of  optimal  orthonormal  (paraunitary)  filter  banks  [1,  2,  3,  4],  This  connection  is  made  for 
the  case  where  the  filters  are  allowed  to  be  ideal.  More  recently  a  number  of  authors  have  considered  the  FIR 
energy  compaction  problem  for  the  two-channel  case  [1,  5,  6,  7,  8,  9,  10]  and  for  the  M-channel  case  [11]. 
The  FIR  compaction  filters  have  applications  in  many  areas  including  data  compression,  signal  analysis,  signal 
modeling,  and  data  transmission  [1,  3,  12],  An  M-channel  FIR  compaction  filter  can  be  considered  as  one  of 
M  filters  of  a  maximally  decimated  M-channel  FIR  orthonormal  filter  bank.  Hence,  one  can  view  the  problem 
as  compaction  of  most  of  the  signal  energy  into  one  channel  of  an  orthonormal  filter  bank. 

In  this  paper  we  consider  some  new  design  techniques  for  optimum  FIR  compaction  filters.  We  give  analytical 
solutions  in  the  two-channel  case  for  a  class  of  random  processes.  Some  properties  of  optimum  FIR  compaction 
filters  and  corresponding  gains  are  also  considered.  Detailed  outline  is  provided  in  Sec.  1.4. 

1.1.  Notations  and  Terminology 

1.  Bold  faced  upper  and  lower  case  letters  represent  matrices  and  vectors  respectively. 

2.  X(z)  and  X(ejuJ)  stand  for  z -transform,  and  Fourier  transform,  respectively  of  a  sequence  x(n).  The 
notation  X(z )  denotes  the  z-transform  of  z*(-n)  where  *  stands  for  complex  conjugation.  If  i(n)  is  real, 
then  X(z )  =  Xiz-1).  Notice  that  X(z)  =  X'(l /z*),  and  the  Fourier  transform  of  x'(-n)  is  X’(e^). 

3.  The  symbols  J.  M  and  t  M  denote  M-fold  decimation  and  expansion  as  defined  in  [13].  The  notation 
X(z)\im  denotes  the  z-transform  of  the  decimated  sequence  x(Mri). 

4.  Nyquist(M)  property.  A  sequence  x(n)  is  said  to  be  Nyquist(M)  if  x{Mn)  =  <5(n)  or  equivalently 

X(z)\im  =  1-  This  can  be  rewritten  in  the  form  [13]: 

M— 1 

£2  X(zWk)  =  M  (!) 

k= 0 

where  W  =  .  In  an  M-channel  orthonormal  filter  bank  {Hk(z)},  each  \Hk{e^)\2  is  Nyquist(M). 

So  the  integer  M  is  often  referred  to  as  the  number  of  channels. 

5.  The  notation  xL(n )  stands  for  a  periodic  sequence  with  periodicity  L.  If  there  is  a  reference  to  a  finite 
sequence  x(n)  as  well,  then  it  is  to  be  understood  that  xL(n)  is  the  periodical  expansion  of  x(n)  with 
period  L,  i.e.,  xL(n)  =  £“-o 0x(n  +  Li).  The  Fourier  series  coefficients  of  xL(n)  is  denoted  by  XL(k). 

6.  For  L  a  multiple  of  M,  a  periodic  sequence  xl{ n)  is  said  to  be  Nyquist(M)  if 

OO 

x l(Mti)  =  Sfcin)  —  ^(n  +  ^ 
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where  K  -  L/M.  The  equivalent  form  of  this  property  in  terms  of  the  Fourier  series  coefficients  XL{k)  is 

M—l 

J2XL(k  +  iK)  =  M,  k  =  0,...,K-l  (3) 

<=o 

(see  Lemma  2  in  Sec.  III). 

7.  Positive  definite  sequences.  Let  a  sequence  {x(n),n  =  0, ...  ,1V}  be  given  and  let  P  be  the  Hermitian 
Toeplitz  matrix  whose  first  row  is  [x(0)  x(l)  ...  x(N)],  The  sequence  {x(n),n  =  0,  ...,1V}  is  called 
positive  (negative)  definite  or  semidefinite  if  P  is  positive  (negative)  definite  or  semidefinite  respectively. 
Let  [a(0)  a(l)  . . .  a(N)]T  denote  the  corresponding  eigenvector.  Then  the  filter  A(z)  =  £n=o  a(n)2_" 
will  be  called  a  maximal  eigenfilter  of  P.  If  we  consider  the  minimum  eigenvalue  instead,  we  will  call 
the  corresponding  filter  a  minimal  eigenfilter  of  P . 


1.2.  The  FIR  energy  compaction  problem 

A  filter  H(z)  of  order  N  will  be  called  a  valid  compaction  filter  for  the  pair  (M,  N)  if  the  product  G(z)  = 
H(z)H(z)  is  Nyquist(M).  We  will  refer  to  G(z)  as  the  product  filter  corresponding  to  H(z).  Conversely, 
G(z)  is  the  product  filter  of  a  valid  compaction  filter  for  the  pair  ( M,N )  if  it  is  of  symmetric  order  N,  that  is 
G(z)  =  2(n)z_n  and  i1;  satisfies  the  following  conditions: 

g{Mn )  =  6(n)  and  G(eju)  >  0.  (4) 


Now  consider  Fig.  1  where  H(z)  is  an  FIR  filter  of  order  N  applied  to  an  input  x(n)  which  is  a  zero-mean 
WSS  random  process  with  the  power  spectral  density  5Ii(eJW).  The  output  of  the  filter  is  decimated  by  M  to 
produce  y{n).  The  optimum  FIR  energy  compaction  problem  is  to  find  a  valid  compaction  filter  H{z)  for  the 
pair  (M,  N)  such  that  the  variance  <r£  of  y{n)  is  maximized.  Since  decimation  of  a  WSS  process  does  not  alter 
its  variance,  we  have 


a2y  =  f  \H(eni2Sxx(en^ 

J  — JT 

We  will  consistently  use  the  notation  G(ejLJ)  =  \H(e^)\2.  We  define  the  compaction  gain  as 

^  comp\lvl  ? iV  /  o  C 

ax  J-n^ xx\eJ  >2^ 


(5) 

(6) 


where  is  the  variance  of  x(n).  The  aim  therefore  is  to  maximize  the  compaction  gain. 

Special  cases.  There  are  two  extreme  cases  worth  examining:  the  case  where  N  <  M  and  the  case  where  ideal 
filters  are  allowed.  In  the  first  case,  the  condition  g(Mn)  =  6(n)  is  the  same  as  g( 0)  =  1.  This  is  equivalent 
to  saying  that  H(eju)  has  unit  energy.  Hence,  by  Rayleigh’s  principle  [14],  the  best  filter  is  the  maximal 
eigenfilter  of  the  Hermitian  Toeplitz  matrix  P  whose  first  row  is  [r(0)  r(l)  . . .  r(N)]  where  r(n)  =  r*(-n) 
is  the  autocorrelation  sequence  of  x(n).  The  corresponding  compaction  gain  is  the  maximum  eigenvalue  of 
P  and  will  be  called  the  KLT  gain.  The  second  case  has  been  studied  in  [4]  where  the  author  describes  a 
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constructive  solution  that  involves  comparison  of  a  set  of  values  of  the  power  spectral  density  at  each  frequency 

u>  and  assignment  of  certain  values  to  H  (eJW)  accordingly. 

In  the  FIR  energy  compaction  problem,  we  do  not  have  the  flexibility  of  assigning  values  to  H(e?u)  inde¬ 
pendently  for  each  w.  This  is  because  H(eju)  is  determined  by  its  IV  +  1  frequency  samples.  For  N  >  M, 
the  problem  is  not  an  eigenfilter  problem  either,  as  the  condition  g{Mn)  =  S(n)  implies  more  than  the  simple 
unit-energy  condition.  In  Sec.  Ill  we  will  introduce  a  suboptimal  method  called  the  window  method  for  design 
of  such  filters.  Interestingly  enough,  this  design  method  involves  two  stages  that  can  be  associated  with  the 
above  two  extreme  cases.  While  the  method  is  not  optimal,  it  produces  compaction  gains  very  close  to  the 
optimum  ones  especially  for  high  filter  orders. 

1.3.  Previous  work 

The  major  motivation  for  studying  compaction  filters  is  their  connection  to  optimal  subband  coding  problem 
which  has  been  widely  studied  [1,  2,  4,  15,  16,  17,  18].  When  the  ideal  filters  are  allowed,  the  optimization  of  a 
maximally  decimated  M-channel  orthonormal  filter  bank  for  a  given  input  statistics  has  been  solved  [2,  3,  4], 
and  the  biorthogonal  case  is  addressed  in  [18].  Optimization  of  FIR  orthonormal  filter  banks  is  analytically  more 
difficult,  and  some  numerical  methods  have  been  developed  for  the  two-channel  case  [1,  5,  9,  10]  and  for  the 
M-channel  case  [11].  Lattice  parametrization  is  utilized  in  [1]  and  the  parameters  are  iteratively  optimized.  An 
optimum  M-channel  FIR  compaction  filter  is  designed  by  linear  programming  in  [11]  and  then  an  orthonormal 
filter  bank  is  constructed  in  some,  optimal  sense  using  the  remaining  degrees  of  freedom.  The  FIR  energy 
compaction  problem  has  been  considered  by  several  authors.  There  are  different  approaches  to  the  problem 

summarized  below. 

1.  Eigenfilter  method.  In  [19],  the  authors  design  one  filter  of  an  M-channel  orthonormal  filter  bank 
using  the  so-called  eigenfilter  method.  The  objective  in  their  design  is  to  have  a  good  frequency  response. 
However,  one  can  modify  the  technique  to  incorporate  the  input  statistics.  This  can  be  done  by  using  the 
psd  Sxx(ejuJ)  as  a  weighting  function  in  the  optimization.  The  paper  also  discusses  how  to  design  a  good 
orthonormal  filter  bank  using  the  remaining  degrees  of  freedom.  In  [11]  the  authors  show  how  to  use  this 
idea  for  the  statistical  optimization  of  orthonormal  filter  banks. 

2.  Linear  programming.  If  one  considers  the  problem  of  finding  the  product  filter  G(z )  corresponding  to 
the  optimum  compaction  filter  H(z),  then  it  can  be  formulated  as  a  linear  programming  problem.  This 
has  been  done  recently  by  Moulin  [9,  11].  For  a  brief  description  of  the  formulation,  the  reader  is  referred 
to  Sec.  2.1.  The  compaction  filter  H(z)  is  obtained  from  G(z)  by  spectral  factorization. 

3.  Quadratic  constrained  optimization  method.  If  one  formulates  the  FIR  energy  compaction  problem 
in  terms  of  the  compaction  filter  impulse  response  h(n),  rather  than  the  product  filter  g(n),  then  it 
becomes  a  quadratic  constrained  optimization  problem  with  a  quadratic  objective.  In  this  method  no 
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spectral  factorization  is  involved.  This  method  has  been  applied  for  the  two-channel  case  by  Caglar 
et  al.  [5],  The  general  M -channel  case  has  been  considered  in  a  completely  different  context  in  the 
communications  literature  by  Chevillat  and  Ungerboeck  [12]  where  an  explicit  algorithm  flow-chart  is 
provided.  In  [20],  a  filter  with  good  frequency  response  is  designed  using  a  similar  approach.  This  is  a 
valid  compaction  filter  according  to  our  definition  (Sec.  1.2),  and  the  extension  to  the  optimal  compaction 
filter  design  is  straightforward. 


4.  Analytical  methods.  Aas  et  al.  [21]  worked  on  a  closely  related  problem  for  the  two-channel  case  even 
though  they  have  not  explicitly  considered  the  energy  compaction  problem.  These  authors  have  optimized 
a  real-coefficient  FIR  filter  H(z)  such  that  H{z)H{z~l)  is  Nyquist(2)  (that  is,  it  is  a  valid  real-coefficient 
compaction  filter),  and 

(7> 

is  maximum.  The  elegance  of  the  method  in  [21]  lies  in  the  fact  that  no  iterative  numerical  optimization  is 
involved.  Based  on  the  fundamentals  of  Gaussian  quadrature,  the  authors  were  able  to  obtain  an  analytical 
method  to  identify  the  unit-circle  zeros  of  H  (z)  which  uniquely  determine  it.  In  our  paper,  this  method  wall 
be  referred  to  as  analytical  method.  In  Sec.  V  we  present  extensions  of  the  analytical  method.  While 


the  original  method  primarily  addresses  conventional  lowpass  filter  design,  we  will  show  how  to  adapt  the 
idea  for  the  case  of  FIR  compaction  filter  design  for  a  given  power  spectrum.  Interestingly  enough,  we 
shall  show  that  the  analytical  method  is  related  to  the  well-known  line-spectral  theory  in  signal  processing 
society  [22].  We  also  mention  here  the  work  of  Usevitch  and  Orchard  [8]  where  an  analytical  expression 
for  the  compaction  gain  is  presented  for  IV  =  3  and  M  =  2. 


1.4.  Summary  of  the  new  results  and  outline  of  the  paper 

1.  Modification  of  linear  programming  method.  In  Sec.  2.2  we  give  a  simple  procedure  to  ensure  that 
linear  programming  solutions  guarantee  the  nonnegativity  requirement  on  G(eJ~').  Since  linear  program¬ 
ming  normally  guarantees  the  nonnegativity  only  at  a  selected  number  of  frequencies,  it  is  important  to 
make  a  modification  to  be  able  to  obtain  a  solution  so  that  G(ejw)  >  0,  Viu. 

2.  Multistage  (IFIR)  extensions.  In  Sec.  2.3  we  present  the  so-called  IFIR  extension  of  the  linear  pro¬ 
gramming  technique  to  design  FIR  compaction  filters  in  two  stages.  That  is,  the  filter  H ( z )  is  represented 
in  the  form  H0(z)Hi(zm°)  where  Mo  is  a  factor  of  M  and  the  much  smaller  filters  Hq(z)  and  H\{z)  are 
optimized.  While  theoretically  suboptimal,  the  system  H0(z)Hi(zM°)  offers  a  much  higher  filter  order  for 
a  fixed  number  of  coefficients.  This  makes  the  design  as  well  as  implementation  very  efficient,  similar  to 
the  case  of  IFIR  filters  in  filter  design  practice  [13,  23].  In  Sec.  2.4  we  consider  a  slight  variation  of  the 
configuration  where  the  Nyquist(M)  property  is  theoretically  assured  and  any  compaction  filter  design 
technique  can  be  used  to  design  individual  filters.  The  multistage  method  is  also  motivated  by  the  fact 
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that  in  the  case  of  ideal  filters  (i.e.,  when  there  is  no  order  constraint),  this  approach  does  not  result  in  a 
loss  of  compaction  gain  as  proved  in  [24]. 

3.  Window  method.  A  new  and  efficient  design  technique  is  introduced  in  Sec.  Ill,  called  the  window 
method  for  the  design  of  FIR  compaction  filters.  For  a  given  input  power  spectrum,  this  method  yields  a 
suboptimal  solution  which  is  very  close  to  the  optimal  solution  especially  for  large  filter  orders.  The  window 
method  has  the  advantage  that  no  optimization  tools  or  iterative  numerical  techniques  are  necessary.  The 
solution  is  generated  in  a  finite  number  of  elementary  steps,  the  crucial  step  being  a  simple  comparison 
operation  on  a  finite  frequency  grid.  Combined  with  the  fact  that  the  solution  is  close  to  optimal,  the 
method  offers  an  attractive  alternative  to  linear  programming.  In  fact,  we  will  show  in  Sec.  3.2  that  there 
is  a  connection  between  the  two  methods. 

4.  Properties  of  optimum  FIR  compaction  filters  and  gains.  In  Sec.  IV  we  examine  some  properties 
of  optimum  FIR  compaction  filters  and  corresponding  maximum  compaction  gains.  For  example,  we  give 
a  necessary  condition  for  a  compaction  filter  to  be  optimum  for  a  given  power  spectrum.  It  is  shown  that 
an  optimum  FIR  compaction  filter  continues  to  be  optimum  if  the  power  spectrum  is  modified  in  certain 
ways.  The  behaviour  of  the  optimum  compaction  gain,  denoted  by  Gopt(M,  N),  as  a  function  of  M  and  N 
for  a  given  power  spectrum  is  investigated.  We  give  some  useful  lower  and  upper  bounds  on  Gopt(M,N). 

5.  Analytical  method.  We  consider  the  FIR  energy  compaction  problem  analytically  for  the  two-channel 
case  (Sec.  V).  This  is  possible  for  a  certain  class  of  random  processes  only.  It  can  be  regarded  as  a 
generalization  of  the  technique  in  [21].  For  the  algorithm  to  be  applicable,  a  certain  sequence  derived 
from  the  odd  autocorrelation  sequence  has  to  be  a  positive  or  negative  definite  sequence.  We  characterize 
classes  of  random  processes  for  which  this  is  the  case  and  therefore  the  method  is  applicable.  As  examples, 
we  give  analytical  solutions  for  MA(1)  and  AR(1)  processes. 

1.5.  Connection  between  energy  compaction  filters  and  optimum  orthonormal  filter  banks 

For  the  two-channel  orthonormal  subband  coder  shown  in  Fig.  2,  the  coding  gain  expression  takes  the  form  [4] 

r  -  (8) 

where  the  second  line  follows  from  the  condition  2 o\  =  +  <  imposed  by  orthonormality.  Here  o*  is 

the  variance  at  the  output  of  H{(z).  Maximizing  the  coding  gain  is  therefore  equivalent  to  maximizing  (or 
minimizing)  cr^Q  under  the  Nyquist(2)  constraint 

H0(z)H0(z)\i2  =  l.  <9) 

Notice  that  the  above  argument  holds  even  if  H0(z)  is  constrained  to  be  a  finite  order  (HR  or  FIR)  transfer 
function.  In  the  M- channel  case,  the  coding  gain  is  still  the  AM/GM  ratio  of  subband  variances,  but  it  cannot 
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be  expressed  in  terms  of  a  single  subband  variance  a2Xg .  It  can  however  be  shown  that  if  the  filter  orders  are 
unconstrained,  then  the  analysis  filters  are  optimal  compaction  filters  for  appropriate  power  spectra  derived 
from  the  input  [4].  For  the  finite  order  case  and  arbitrary  M,  optimal  compaction  filters  are  still  of  interest 
beacuse  of  the  large  coding  gain  obtainable  from  them. 

II.  LINEAR  PROGRAMMING  METHOD  AND  IFIR  DESIGNS 

The  use  of  linear  programming  method  in  compaction  filter  design  was  proposed  in  [9],  and  is  reviewed  in 
Sec.  2.1.  The  technique  yields  a  solution  G(e>“)  that  is  Nyquist(M)  and  is  nonnegative  on  a  certain  discrete 
grid  of  specified  frequencies.  Note  that  after  finding  G(z),  one  needs  to  spectrally  factorize  it  to  find  the 
compaction  filter  H(z).  This  step  will  succeed  only  if  G(e^)  >  0,  Vw.  In  Sec.  2.2,  we  propose  a  simple 
procedure  to  guarantee  the  nonnegativity  of  the  resulting  solution  for  all  frequencies.  We  then  give  extensions 
of  the  technique  for  the  case  of  multistage  (IFIR)  compaction  filters  in  Sec.  2.3.  We  will  also  consider  a  special 
IFIR  configuration  in  Sec.  2.4  which  has  certain  advantages. 

2. 1.  Review  of  the  linear  programming  method 

Assume  that  the  input  process  x(n )  is  real.  The  output  variance  can  be  written  as 

&y  =  r(0)  +  2  g{n)r(n)  (10) 

n=l 

Let  g d  and  rd  be  the  vectors  formed  by  the  nonzero  components  of  g(n)  and  r(n)  for  n  =  1, ...  IV.  That  is, 

g d  =  [g(l)  g( 2)  . . .  g(M  -  1)  g(M  +1)  ...  g(N)]T,  rd  =  [r(l)  r(2)  . . .  r(M  -  1)  r(M  +  1)  ...  r(N)f  (11) 

Then  (10)  can  be  written  as  oj  =  r(0)  +2rjg<i.  This  incorporates  the  Nyquist (M)  condition  but  not  the 
nonnegativity  constraint  in  (4).  Let  Cd(u)  =  [cos(w)  cos(2w)  ...  cos((M  —  l)ui)  cos((Af  +  l)w)  ...  cos(IVu>)] 
Then  G(ejtJ)  =  1  +  2cd{uj)gd.  Hence  the  problem  is  equivalent  to  the  following: 

maximize  rdgd  subject  to  cd(w)gd  >  —0.5,  Vw  €  [0,tt]  (12) 

This  type  of  problem  is  typically  classified  as  semiinfinite  linear  programming  [9].  By  discretizing  the  frequency, 
one  reduces  this  to  a  well  known  standard  linear  programming  problem.  This  discretization  however  needs  to 
be  done  with  care.  The  resulting  G(ej“)  can  go  negative  in  between  the  discrete  frequencies.  To  avoid  this,  one 
can  put  some  tolerance  in  the  inequality  (12).  Our  experience  however  was  that  even  with  a  high  tolerance, 
the  resulting  G(z)  had  single  (rather  than  double)  unit-circle  zeros  and  G(eJ“)  was  not  nonnegative.  One  way 
to  overcome  the  difficulty  is  to  numerically  determine  the  zeros  of  G(e]U)  and  to  merge  the  pairs  of  zeros  that 
are  very  close  to  the  unit-circle  into  double  unit-circle  zeros.  This  requires  determining  the  roots  of  G(z)  that 
are  in  the  vicinity  of  the  unit-circle.  This  can  be  done  by  looking  at  the  Fourier  transform  G(eJU)  using  a 
large  number  of  frequency  points.  Yet  another  way  is  to  “lift”  G(eju)  by  increasing  5(0)  relative  to  the  other 
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coefficients  (since  g(0)  has  to  be  1,  in  effect  we  scale  g(n)  for  n  #  0  by  a  constant  that  is  slightly  less  than  1). 
In  the  next  section  we  propose  a  third  technique  to  overcome  the  difficulty  without  having  to  locate  any  zeros 

or  the  minimum  of  G(eJU’). 

2.2.  Windowing  of  the  linear  programming  solution 

Consider  the  periodical  expansion  gL(n)  of  the  linear  programming  solution  where  L  is  the  number  of 
discrete  uniform  frequencies  {uk}  used  in  the  design  process.  Assume  that  L  >  2 N.  Linear  programming 
assures  that  G(eju)  is  nonnegative  at  the  frequencies  {uk}.  Hence  the  Fourier  series  coefficients  GL(k)  of  gL(n) 
are  nonnegative.  Now  if  we  consider  the  product 

w(n)gL(n)  (13) 

where  w(n)  is  a  symmetric  window  of  order  K  <  L  -  N  (length  2 K  +  1),  the  resulting  Fourier  transform 
is  nonnegative  provided  w(n)  has  nonnegative  Fourier  transform  W(ej“).  This  is  depicted  in  Fig.  3.  The 
reason  follows  from  the  fact  that  the  Fourier  transform  of  w(n)gL(n)  is  a  weighted  sum  of  shifted  versions  of 
W(ejw)  with  nonnegative  weights.  For  maximum  compaction  gain,  the  symmetric  order  of  w(n)  is  chosen  to 
be  maximum,  namely  K  =  L  -  N  -  1.  Note  that  when  L  =  2N,  we  have  gL(N)  =  2 g(N).  One  can  use  a  fixed 
window  like  a  triangular  window  as  depicted  in  the  figure  and  get  a  satisfactory  compaction  gain.  However 
one  can  always  optimize  the  window.  The  optimization  of  the  window  given  the  periodic  sequence  gL(n)  is 
discussed  in  Sec.  Ill  where  we  show  that  the  optimum  w(n)  is  the  product  filter  of  the  maximal  eigenfilter  of 
the  K  x  K  Hermitian  Toeplitz  matrix  formed  by  the  product  r(n)gL(n).  Since  the  symmetric  window  order  K 
is  very  high  in  linear  programming  designs,  we  suggest  to  use  a  triangular  window  rather  than  optimizing  the 
window.  The  performance  loss  is  negligibly  small. 

Example  1.  Let  the  input  be  psd  be  as  in  Fig.  4  and  let  N  =  65  and  M  =  2.  In  the  same  figure,  we  plot  the 
magnitude  square  \H(eju)\2  of  the  compaction  filter  H(z)  designed  by  the  linear  programming  method.  The 
number  of  frequencies  used  in  the  design  process  was  L  =  512.  We  have  used  triangular  window  of  symmetric 
order  K  =  L  —  N  -  1  =  446  and  found  that  the  resulting  compaction  gain  is  GCOmp( 2,65)  =  1.8698.  If  we 
optimize  the  window  the  compaction  gain  becomes  1.8744.  One  can  verify  that  the  compaction  gain  of  the  ideal 
(infinite  order)  compaction  filter  is  1.8754. 

2.3.  IFIR  designs  using  linear  programming 

In  the  design  of  narrowband  FIR  lowpass  filters  for  conventional  applications,  it  is  possible  to  decompose  the 
transfer  function  H(z)  into  the  form  H(z)  =  where  H0(z)  and  H,{z)  have  significantly  smaller 

lengths  than  H(z).  While  H(z)  is  theoretically  suboptimal  compared  to,  e.g.,  an  equiripple  solution,  it  has 
the  advantage  of  actually  requiring  fewer  active  multiplier  elements  (because  of  the  zero-valued  coefficients  in 
Hi(zm°)).  This  technique,  called  the  IFIR  technique  [23]  has  also  been  extended  to  bandpass  and  multiband 
filters  in  the  past,  and  is  in  fact  related  to  multistage  design  of  interpolation  filters  [25].  The  method  offers 
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significant  economy  both  in  terms  of  design  time  and  implementation  complexity.  For  the  case  of  compaction 
filter  design,  a  similar  decomposition  proves  to  be  valuable  as  we  shall  now  demonstrate. 

Assume  M  =  M0Mi.  Consider  Fig.  5(a)  where  H0(z)  and  Hi(z)  are  such  that  the  equivalent  filter  in 
Fig.  5(b)  H(z)  =  H0(z)Hi(zMo)  is  a  valid  compaction  filter  for  the  pair  (M,N).  The  problem  is  to  optimize 
the  pair  of  filters  H0(z)  and  Hi(z)  for  maximum  compaction  gain.  If  we  fix  one  of  the  filters,  it  will  be  shown 
that  the  design  of  the  other  can  be  formulated  as  a  linear  programming  problem.  We  will  describe  the  details 
of  how  to  find  Hx{z)  for  a  fixed  H0(z)  and  vice  versa,  in  an  iterative  manner. 

Let  G0(z)  =  H0{z)H0(z-1),  Gx(z)  =  Hi(z)Hi{z~%  and  G(z)  =  H(z)H(z~l)  with  impulse  responses 
ffo(n),  91  (n),  and  g(n)  respectively.  Denote  the  orders  of  H0(z),  Hi (z),  and  H(z)  by  N0,  Nlt  and  N  respectively. 
Note  that  N  =  M0Ni  +  N0-  Define 

go  =  bo(0)  go(l)  •  •  •  9o(N0))T,  gi  =  bi(0)  9i0)  •  •  •  S  =  fo(°)  9<X)  ■  ■  ■  9Wf-  (14) 

Optimization  of  H\{z)  for  a  given  Hq{z)‘.  We  have  G(z)  =  G0(z)Gi(zm°).  LetG0  bethe  (27V+l)x(2M07V1  + 
1)  convolution  matrix  formed  by  g0{n).  Taking  into  account  the  symmetries  and  the  fact  that  Gi(zm°)  has 
nonzero  components  only  for  multiples  of  M0,  we  can  write  g  =  A0gi,  where  A0  is  an  (TV  +  1)  x  (Nx  +  1) 
matrix  that  is  obtained  from  G0.  Now,  the  Nyquist(M)  constraint  requires  that  if  we  decimate  g  by  M 
we  should  get  e0  =  [1  0  ...  0]T.  Let  B0  denote  the  matrix  that  is  obtained  by  taking  every  Mth  row 
of  A0.  Then  we  should  have  B0gi  =  e0.  To  force  the  nonnegativity  constraint  on  Gi(e,w),  let  c0(w)  = 
[1  2  cos(w)  2  cos(2w)  . . .  2cos(W1w)]T.  Then  the  constraint  Gi(e^)  >  0  becomes  cj (w)gl  >  0,  Vw  £  [0,tt].  If 
r  _  ^r(o)  2r(l)  . . .  2r(iV)]T,  the  objective  is  to  maximize  rTg  =  rTA0gi.  Hence  we  have  reduced  the  problem 

to  the  following: 

maximize  rggi,  subject  to  B0gi  =  e0,  and  cj (w)gi  >  0,  Wui  £  [0,7r]  (15) 

where  r0  =  Aq  r.  Hence  a  standard  linear  programming  algorithm  can  be  applied  once  a  set  of  freqencies  is 
chosen  for  the  inequality  constraint. 

Optimization  of  H0(z)  for  a  given  i?i(z):  Similarly,  one  can  reduce  the  problem  of  finding  the  best  H0(z) 
for  a  given  Hi(z)  to  the  following  linear  programming  problem: 

maximize  rf  g0,  subject  to  Bigo=e0,  and  cf  (oj)g0  >  0,  Vw  £  [0, 7r]  (16) 

where  rx  =  Afr,  ci(u >)  =  [12 cos(u>)  2 cos(2w)  ...  2 cos (N0uj)]t .  The  (N  + 1)  x  (N0 + 1)  matrix  Ax  is  obtained 
from  the  (2N  +  1)  x  (2N0  +  1)  convolution  matrix  formed  by  51  (n)  by  taking  the  symmetries  into  account  and 

the  matrix  Bi  is  obtained  by  taking  every  Mth  row  of  Ai. 

One  can  iterate  between  the  above  two  optimization  steps  until  there  is  no  significant  change  in  the  com¬ 
paction  gain.  The  initial  choice  of  g0(n)  can  significantly  affect  the  resulting  compaction  gain.  According  to  our 
design  experience  if  go(n )  is  chosen  to  be  a  triangular  sequence,  the  compaction  gam  at  the  end  of  the  iteration  is 
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very  good.  The  filters  g0(n)  and  gi(n)  which  result  from  the  iteration  should  spectrally  be  factorized  to  identify 
Ho(z)  and  Hx(z).  This  step  will  be  successful  only  if  the  solutions  are  such  that  Go(eJU)  >  0  and  Gi(ej  )  >  0 
for  all  u.  If  this  is  not  the  case,  we  can  force  it  by  use  of  windowing  on  g0(n)  and  51 0)  35  described  in  Sec.  2.2. 
If  this  is  done  then  the  product  filter  G0(z)Gl(zM<>)  will  not  be  exactly  Nyquist(M).  In  the  next  subsection  we 
show  how  to  overcome  this  problem. 

Example  2.  Let  us  design  IFIR  compaction  filters  for  the  pair  (M,  N)  =  (36,65),  and  for  the  input  process 
whose  psd  is  given  in  Fig.  4.  Let  M0  =  9  and  M,  =  4,  and  let  N0  =  11  so  that  N,  =  6.  The  number  of  frequencies 
used  in  the  designs  is  L  =  1024.  Starting  with  a  triangular  sequence  for  g0(n),  the  algorithm  converges  m  a 
few  steps.  We  windowed  the  resulting  solutions  g0(n)  and  9l(n)  with  triangular  windows  of  symmetric  orders 
L-No-l  and  L-Ni-1  respectively.  The  final  product  filter  was  not  exacly  Nyquist(M)  because  it  was  found 
that  g{ 36)  ~  -0.0018  ^  0.  The  final  compaction  gain  was  5.1444.  If  we  design  a  compaction  filter  of  order  18 
directly  (i.e.,  not  using  IFIR  technique),  the  compaction  gain  is  4.4225.  This  corresponds  to  a  compaction  filter 
with  the  same  number  of  active  multipliers,  namely  19.  If  we  design  a  compaction  filter  of  order  65  directly  (66 
active  multipliers),  then  the  resulting  compaction  gain  is  7.2337. 

2.4.  A  Particular  IFIR  configuration 

In  Fig.  5,  if  G0(z)  is  Nyquist(M0)  and  Gi(z)  is  Nyquist(M!),  it  can  be  verified  that  G(z)  given  by 
Gq(z)Gi(zMo)  is  Nyquist(M).  Now,  let  us  fix  H0(z)  to  be  a  valid  compaction  filter  for  the  pair  (N0,M0). 
Referring  to  Fig.  6(a),  the  best  Hi(z)  is  the  optimum  compaction  filter  for  (NuMi),  and  for  the  input  x0(n) 
which  has  the  psd  5I0l0(z)  =  (g0(*)S**(z))|w  Similarly,  if  Hy(z)  is  a  fixed  compaction  filter  for  the  pair 
then  we  can  redraw  the  configuration  as  in  Fig.  6(b)  and  therefore  the  best  H0(z)  is  the  optimum 
compaction  filter  for  ( N0,M0 ),  and  for  the  input  *i(»)  which  has  the  psd  SXlXl(z)  =  Gi(zM°)Sxx(z).  One  can 
design  the  compaction  filters  H0(z )  and  Hx(z)  iteratively  using  any  of  the  known  techniques.  Hence,  one  can 
use  the  linear  programming  technique  as  well  as  any  other  technique  like  the  window  method  to  be  introduced 
in  Sec.  III.  Also  note  that  if  the  ideal  filters  are  allowed,  this  multistage  configuration  has  no  loss  of  generality 

as  shown  in  [24]. 

Example  3.  Let  the  setup  be  the  same  as  in  Example  2.  We  have  designed  the  compaction  filters  H0(z ) 
and  Hi(z)  iteratively  using  the  standard  linear  programming  procedure  as  in  Example  1.  We  have  started 
with  Hx(z)  =  1.  The  first  compaction  filter  H0(z)  is  therefore  the  optimal  compaction  filter  for  the  pair 
(Mq,  No)  =  (9,11)  for  the  original  autocorrelation  sequence.  We  have  windowed  the  final  product  filters  as 
we  did  in  Example  1  to  guarantee  the  nonnegativity.  The  resulting  overall  compaction  gain  is  4.9432.  This  is 
slightly  smaller  than  the  overall  compaction  gain  5.1444  in  Example  1.  However,  the  resulting  overall  filter  here 
is  exactly  Nyquist(Af )  unlike  the  case  of  Example  2. 

III.  WINDOW  METHOD 
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In  this  section  we  will  describe  a  new  method  to  design  compaction  filters  for  general  M.  The  input  process 
might  be  real  or  complex  and  its  psd  may  have  any  shape.  The  technique  is  quite  simple  while  the  resulting 
compaction  gains  are  very  close  to  the  optimum  ones  especially  for  relatively  high  filter  orders. 

3.1.  Derivation  of  the  window  method 

The  idea  behind  the  method  is  to  represent  the  impulse  response  of  the  product  filter  G{z)  in  the  form 

g(n)  =  w(n)fL(n),  (17) 

where  w(n )  and  h(n)  are  conjugate  symmetric  (i.e.,  w(n)  =  w*(-n),  fL(n)  =  fi(~n))>  411(1  w(°)  ~  M0)  =  1 
(see  Fig.  7).  The  window  w(n)  has  the  same  length  as  g(n),  namely  2N  4-  1  and  the  sequence  /i,(n)  is  periodic 
with  period  L  =  KM  >  2 N  for  some  K.  Evidently,  only  one  period  of  fL{n)  matters  in  (17).  Define  the  Fourier 
series  coefficients  of  /t(n)  as 

FL(k)  =  £  h(n)WkLn,  WL  =  e-»«'L  (18) 

n=0 

This  is  periodic  with  the  same  period  L  (the  first  period  (Ft(fc),  k  =  0, . . . ,  L  —  1}  being  just  the  DFT  of  the 
sequence  {/z,(n),n  =  0, . . . ,  L  -  1}).  We  have  the  following  observation: 

Lemma  1.  Consider  the  representation  (17)  for  g(ri)  where  w(n)  and  } iXn)  are  as  explained  above.  If  the 
Fourier  transform  W{eiu)  of  w(n)  is  nonnegative  for  all  w,  if  the  Fourier  series  coefficients  FL(k)  of  }l{ n)  are 
nonnegative  for  all  k,  and  if  /i,(n)  is  Nyquist(M)  then  G(z )  is  product  filter  of  a  valid  compaction  filter  for  the 
pair  (M,  N).  That  is,  g(Mn)  =  S(n)  and  G(eju)  >  0. 

Proof.  It  is  readily  verified  that  G{e?u)  =  j  XXk--o  FL(k)W(e^ul  L  ^)-  Since  Fi(k)  >  0  and  W(eJ  )  >  0,  it 
follows  that  G(eju)  >  0.  If  fL(n)  is  Nyquist(M)  then  so  is  g{n)  because  L  >  2N.  ■ 

The  choice  of  L  will  be  discussed  later.  Assume  the  conditions  of  the  lemma  hold  so  that  G(z)  is  product  filter 
of  a  valid  compaction  filter.  If  w(n)  is  fixed,  what  is  the  best  fL{n)  that  maximizes  the  compaction  gain?  To 
answer  the  question  we  first  note  the  following: 

Lemma  2.  A  periodic  sequence  fL(n)  with  period  L  =  KM  is  Nyquist(M),  that  is,  fL(Mn )  =  SK(n),  if  and 

only  if  its  Fourier  series  coefficients  Fi(k)  satisfy  the  following: 

M- 1 

Y,FL{k  +  iK)  =  M,  k  =  0,...,K-l  (19) 

1=0 

A  proof  of  this  can  be  found  in  Appendix  A. 

To  obtain  the  best  fL(n)  let  us  write  the  objective  (5)  in  terms  of  w{n)  and  fL(n): 

N  N 

53  9(n)r*(n)=  w(n)h(n)r*(n) 

n=— iV  n=-N 

Let  f(n)  =  w* (n)r(n)  and  let  SL(k)  be  the  Fourier  series  coefficients  of  its  periodic  expansion  rL(n).  For 
simplicity  assume  that  L  >  2N .  Then  the  objective  can  be  written  as 

o-l  =  £  fL{n)fl(n)  =  \  £  FL(k)SL(k)  (21) 

n=0  k= 0 
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Notice  that  both  FL(k)  and  SL(k)  are  real.  Now  to  incorporate  the  Nyquist(M)  constraint  we  write  the 

preceding  as  K1 

—  +  iK)§L(k  +  iK) 


(22) 


Jfc=0  1=0 


For  a  fixed  k,  let  SL(k  +  i0K)  be  the  maximum  of  the  set  {SL(k  +  iK),  i  -  0, . . . , M  -  1}.  Then  by  (19),  and 
noting  that  FL(k)  >  0,  the  objective  (22)  is  maximized  if  we  assign 

FL(k  +  i0K)  =  M,  and  FL(k  +  UK)  =  0,  l  =  -  1.  (23) 


The  procedure  is  illustrated  in  Fig.  8.  By  repeating  the  process  for  each  k  =  0, . . . ,  K  -  1,  the  Fourier  senes 
coefficients  of  the  best  fL(n)  is  determined.  The  sequence  fL(n)  can  now  be  calculated  by  the  inverse  relation: 

(24) 

^  k= o 

Summary  of  the  window  algorithm 

Assume  a  window  w(n)  of  the  same  symmetric  order  as  g{n)  with  nonnegative  Fourier  transform  has  been 
chosen.  Let  L  =  KM  >  2 N.  Then  the  algorithm  steps  are 


1.  Calculate  SL(k),  the  DFT  coefficients  of  rL(n),  n  =  0, . . . ,  L  -  1,  where  fL(n)  is  the  periodical  expansion 
of  r(n)  =  w*(n)r(n). 

2.  For  each  k  =  0,...,I<  -  1,  determine  the  index  i0  for  which  Sz,(k  +  i0K)  is  maximum,  and  assign 

FL(k  +  i0K)  =  M  and  FL(k  +  itK )  =  0,  l  =  1, . . . ,  M  -  1. 

3.  Calculate  /x,(n)  by  the  inverse  relation  (24).  We  need  only  to  determine  //,(n)  for  n  =  1, . . . ,  N. 

4.  Form  the  product  filter  impulse  response  g(n)  =  w(n)/i(n)  and  spectrally  factorize  G(z)  to  find  the 

compaction  filter  H(z). 

Real  case.  Before  proceeding  to  a  design  example,  consider  the  case  of  real  inputs.  In  this  case,  the  above 
algorithm  can  be  modified  to  produce  real-coefficient  compaction  filters.  Let  us  consider  the  set  of  values 

{SL(k  -I-  iK),  i  =  0, . . . ,  M  -  1}  (25) 


for  each  k  =  0 . K  -  1.  Since  SL(k)  =  Sl(L  -  k)  if  the  process  is  real,  the  above  set  is  equivalent  to 

{Sl(L  —  k  —  iK)  =  Sl(K  —  k  +  K(M  —  1  —  i)),  i  =  0, . . .  ,M  —  1}  (26) 

Hence  in  the  comparison,  we  need  to  consider  only  k  =  0, . . . ,  P  where  P  =  -j  if  #  is  even>  anc*  P  ~  2  lf 

it  is  odd.  Let  SL{k  +  ioK)  be  the  maximum  of  the  set  (25)  for  each  k  =  0, . . . ,  P.  Because  of  the  symmetry 
requirement,  we  need  to  be  careful  in  the  assignments.  There  are  two  cases  to  consider: 

1.  The  index  L-k-  ioK  is  among  the  set  {A:  -I-  iK,  i  =  0, . . . ,  M  -  1}, 
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2.  It  is  not. 


The  first  case  happens  if  and  only  if  2k  mod  K  =  0.  This  happens  if  k  —  0  or  k  —  2  .  We  assign  Fi(k  -f  i$K) 
Fl(L  -k-  i0K)  =  f  if  k  +  i0K  t  f ,  and  FL{k  +  i0K)  =  M  if  fc  +  i0K  =  f .  In  the  second  case,  we  assign 
F^Jfc  +  i0F)  =  Fl(L  -k-  i0K)  =  M  ifk  +  i0K  ±  0  and  FL(A:  +  i0K)  =  M\tk  +  i0K  =  0.  In  either  case,  we 
set  the  remaining  values  in  the  set  {Fi(k  +  i M ) }  to  zeros.  This  will  maximize  the  objective  (22),  and  fi (n) 
calculated  by  the  inverse  relation  (24)  is  the  best  sequence  and  it  is  real. 

Summary  of  the  window  algorithm  for  the  real  case 

Assume  a  real  symmetric  window  w(n)  of  order  N ,  with  nonnegative  Fourier  transform  is  given.  Let 
L  =  KM  >  2N  as  before.  Let  P  be  as  explained  above.  Then  the  algorithm  for  the  real  processes  has  the 
following  steps: 

1.  Calculate  SL(k ),  the  DFT  coefficients  of  rL(n),n  =  0, . . . ,  L  -  1,  where  rL{n)  is  the  periodical  expansion 
of  f(n)  =  w(n)r(n). 

2.  For  each  k  =  0, . . .  ,P,  determine  the  index  z0  for  which  SL{k  +  i0K)  is  maximum, 

3.  If  k  +  i0K  =  0  or  k  +  i0K  =  %  then  set  FL(k  +  i0K)  =  M,  else  if  A;  =  0  or  k  =  f ,  then  set  F(k  +  i0K)  = 
F(L  -  Jfc  -  i0K)  =  4f ,  else,  set  F(k  +  i0K)  =  F(L-k-  i0K)  =  M.  Set  all  the  remaining  values  to  zeros. 

4.  Calculate  /L(n)  by  the  inverse  relation  (24).  We  need  only  to  determine  /r(n)  for  n  =  1 

5.  Form  the  product  filter  g{n)  —  w(n)  f l(ji)  and  spectrally  factorize  G (z)  to  find  the  compaction  filter  H (z). 

Optimization  of  the  window.  The  algorithm  produces  very  good  compaction  gains  especially  when  the  filter 
order  is  high  as  we  shall  demonstrate  shortly.  However,  one  can  get  better  compaction  gains  by  optimizing  the 
window  w(n).  Consider  the  representation  (17)  again  and  let  w(n)  and  f l{11)  satisfy  the  conditions  of  Lemma 
1.  If  we  fix  /z,(n),  what  is  the  best  window  w(n)l  The  objective  (5)  can  be  written  as 

al  =  £  Sxx(ejuW{en^  (27) 

where  W(eju)  is  the  Fourier  transform  of  w(n)  and  Sxx(eju)  is  the  Fourier  transform  of  /*(n)r(n)  where 
f(n)  is  one  period  of  f l (n)  centered  at  n  =  0.  Let  W(e^)  =  |.4(e-,^')|2,  where  A(z)  =  )Cn=o  ain)z 
the  spectral  factor  of  W{ejw).  The  only  constraint  on  A{eju)  is  that  it  has  to  have  unit  energy  in  view  of 
w(0)  =  f*K \A(eju)\2$%  =  1.  Hence,  by  Rayleigh’s  principle  [14],  (27)  is  maximized  if  A{z)  is  the  maximal 
eigenfilter  of  P.  The  corresponding  compaction  gain  is  the  maximum  eigenvalue  of  P. 

We  have  described  how  to  optimize  w(n)  given  and  vice  versa.  It  is  reasonable  to  expect  that  one  can 

iterate  and  obtain  better  compaction  gains  at  each  stage.  We  have  observed  that  two  stages  of  iterations  were 
sufficient  to  get  near-optimal  compaction  gains.  We  started  with  a  triangular  window  and  found  that  / i(n)  did 
not  change  after  the  reoptimization  of  the  window.  Notice  that,  the  use  of  an  initial  window  is  not  necessary  if 
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one  is  willing  to  use  a  window  after  finding  fL(n).  However,  in  most  of  the  design  examples  we  considered,  we 
have  observed  that  using  an  initial  window  with  nonnegative  Fourier  transform  (in  particular,  the  triangular 
window)  and  then  reoptimizing  the  window  resulted  in  better  compaction  gains. 

Example  4:  MA(1)  process.  Let  N  =  5,  M  =  4,  r(0)  =  1,  r(l)  =  p,  and  r(n)  =  0,  n  >  1.  Assume  the 
process  is  real  so  that  r(— n)  =  r(n).  Let  the  window  be  triangular,  i.e., 


The  Fourier  transform  of  f(n)  =  w 
in  step  1  are 


n  =  0,  ±1,  - . . ,  ±5 
elsewhere. 


(28) 


(n)r(n)  is  S(eju})  =  1  +  §pcosu>.  Hence,  the  DFT  coefficients  SL(k)  of  f(n) 


Pi  O'jr 

Si,{k)  =  1  + —p  cos— k,  k  —  0, . . . ,  L  —  1. 


(29) 


Now,  assume  L  =  12  >  10,  so  that  K  =  3  and  P  =  1.  So  we  have  the  following  sets  to  consider  in  step  2: 


{SL(  0),  SL{  3),  Si(6),5i(9)},  {SL{  1),  5^(4),  SL(  7),  5l(10)} 


(30) 


5%/3 


which  are  evaluated  below  respectively: 

5  5  5\/3  ,  o  ,  ovo  !  .  /oi\ 

(1  +  3P,1,1~  3P’1^’  ^1  +  -6_p’1-  6Pj1  6  p,1+ ^  ' 

First  assume  p  >  0.  The  maximum  of  the  first  set  is  5jr,(0)  and  the  maximum  of  the  second  set  is  SL(  1).  Hence 

applying  step  3  of  the  algorithm  we  have 

{FL(k),  fc  =  0, . . . ,  L  —  1}  =  {4, 4, 0,0, 0,0, 0,0, 0,0, 0,4}  (32) 


By  the  inverse  relation  (24)  we  calculate  in  step  4: 

{}L{n),n  =  0,..., N}  -  {1, 


1  +  ^  2  1  n  l-y/3, 
3  ’3’3’U’  3  * 


Hence  the  product  filter  g{fi)  —  win) f i{n)  has  been  found,  and 


G(z) 


+  +  ^  + .  +  + 1*-’  +  ^  ^ 


(33) 


(34) 


The  corresponding  compaction  gain  is  1  +  5-^-p  a  1  +  1.5178 p.  An  optimum  compaction  filter  H(z)  is 
obtained  by  spectrally  factorizing  G(z).  Next  consider  the  case  p  <  0.  Let  fL(n)  be  the  corresponding  solution 
with  the  Fourier  series  coefficients  FL(k).  Referring  to  (31),  SL( 6)  in  the  first  set  and  SL(7)  in  the  second  set 
is  maximum.  Hence, 

{FL(k),  k  =  0, . . . ,  L  —  1}  =  (0,0, 0,0, 0,4, 4, 4,0, 0,0,0}  (35) 

which  is  equal  to  FL(k  -  6)  where  FL(k)  is  the  previous  solution.  Hence  fL(n)  =  (-1  )nh(n)  and  therefore 
G(z)  =  G(-z ),  with  the  corresponding  compaction  gain  1  -  An  optimal  compaction  filter  is  H(-z), 

where  H(z)  is  a  solution  for  the  previous  case. 
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For  comparison,  we  have  also  designed  an  optimum  compaction  filter  using  the  linear  programming  technique. 
The  corresponding  compaction  gain  is  approximately  1  +  1.6657|p|.  This  is  achieved  by  using  L  —  512  and  a 
triangular  window  of  symmetric  order  L  —  N  —  1.  The  compaction  gain  of  the  window  method  is  only  slightly 
lower.  Let  us  find  the  improvement  we  can  get  by  optimizing  the  window  when  we  fix  fL(n).  The  compaction 
gain  is  the  maximum  eigenvalue  of  the  6x6  symmetric  Toeplitz  matrix  with  the  first  row  [1  /l(1)  p  0  0  0  0]. 
This  eigenvalue  is  1  +  1.8019/L(l)|p|.  Using  fL(l)  given  in  (33),  the  improved  compaction  gain  is  1  +  1.6410|p| 
which  is  very  close  to  the  linear  programming  compaction  gain  1  +  1.6657|p|. 

Can  we  improve  the  compaction  gain  further  given  this  optimal  window  by  reoptimizing  f i  (n)7  In  this  and 
all  the  other  design  examples  we  considered,  we  used  the  triangular  window  and  then  found  the  optimum  /i(n), 
and  then  reoptimized  w(n)  for  //,  (n).  Interestingly  enough,  the  reoptimization  of  f i  (n)  did  not  change  it! 
Choice  of  the  periodicity  L 

Increasing  L  does  not  necessarily  increase  the  resulting  compaction  gain.  For  example  using  L  —  oo  which 
corresponds  to  using  optimum  ideal  filter  /i,(n)  for  the  autocorrelation  sequence  fit  (n)  does  not  result  in  the 
best  achievable  compaction  gain  using  the  algorithm.  This  is  true  even  if  no  initial  window  w(n )  is  used.  For 
the  above  example,  we  increased  L  to  16  and  found  that  the  compaction  gain  decreased!  When  we  used  the 
ideal  filter  for  /l(u)  which  corresponds  to  L  =  oo,  the  compaction  gain  was  better  than  that  of  the  case  L-  16 
but  worse  than  that  of  the  case  L  =  12. 

The  best  period  L.  Until  this  point  we  assumed  that  L  >  2 N.  With  this  choice  fL{n)  is  the  periodical 
expansion  of  r(n)  with  no  aliasing  (see  Fig.  9(a)).  The  first  period  of  tl(p)  is 


(rx,(n),  n  =  0, . . .  ,L  —  1}  =  (f(0),r(l), . . .  ,t(N), 0, . . . , 0,r*(7V), . . .  ,r*(l)}. 


(36) 


However,  for  the  algorithm  to  successfully  find  a  valid  compaction  filter,  it  is  only  necessary  to  use  a  period  L 
which  is  a  multiple  of  M  and  is  greater  than  N.  This  will  ensure  that  fL(n)  will  be  Nyquist(M).  The  smallest 
such  period  is  M\N/M].  This  choice  however  leads  to  an  additional  symmetry  in  fL(n)  and  according  to  our 
experience,  the  corresponding  compaction  gains  are  not  good.  If  we  use  a  period  L  that  is  the  smallest  multiple 
of  M  such  that  L  >  2 N,  then  we  obtain  very  good  compaction  gains.  This  choice  can  be  compactly  written  as 

L  =  M\2N/M ]  (37) 


If  L  =  2 N,  the  sequence  fii,(n)  has  the  following  first  period: 


{fL(n),  n  =  0, . . . ,  L  -  1}  =  (r(0),  f  (1), . . . ,  f{N)  +  r*(N), . . . ,  f  *  (1)}. 


(38) 


This  is  illustrated  in  Fig.  9(b)  for  a  real  process.  In  this  case,  we  have  fL(N)  =  2 f(N).  This  will  always  be  the 
case  if  M  =  2,  since  L  =  27V  is  a  multiple  of  M. 


3.2.  Connection  between  the  linear  programming  and  window  methods 
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In  both  the  linear  programming  and  window  methods,  we  use  windows  to  assure  the  nonnegativity  of 
G{e?u).  Consider  the  equations  (13)  and  (17).  When  L  is  a  multiple  of  M,  a  periodic  sequence  gL(n)  in  the 
linear  programming  method,  and  a  periodic  sequence  /t(n)  in  the  window  method  are  found  such  that  they 
are  Nyquist (M)  and  their  Fourier  series  coefficients  are  all  nonnegative.  For  L  >  2 N,  two  problems  are  not  the 
same  because  gL(n)  is  necessarily  zero  for  some  n,  while  fL{ri)  can  be  nonzero  for  all  n  (except  n  =  kM,  of 
course).  If  however  L  =  2 N,  then  the  two  problems  are  exactly  the  same!  If  windowing  is  done  in  the  same  way 
in  both  methods,  then  we  see  that  the  resulting  compaction  gains  should  be  the  same.  Hence,  one  can  view  the 
window  method  as  an  efficient  and  noniterative  technique  to  solve  a  linear  programming  problem  when  L  =  2 N. 
If  L  is  inreased,  we  saw  that  the  window  method  does  not  necessarily  yield  better  gams  whereas  this  is  the  case 
for  the  Unear  programming  method  provided  the  window  order  is  increased  as  well.  However,  optimization  of 
the  window  becomes  costly  as  the  order  increases.  If  one  uses  a  fixed  triangular  window  (with  a  high  order)  in 
the  linear  programming,  and  if  the  windows  are  optimized  in  the  window  method,  then  window  method  is  very 
close  and  sometimes  superior  to  the  linear  programming  method  as  we  demonstrate  in  the  following  example. 
Example  5:  Comparison  of  Unear  programming  and  window  methods.  Let  the  input  power  spectrum 
be  as  in  Fig.  4.  In  Fig.  10(a)  we  plot  for  a  fixed  number  of  channels  of  M  =  2,  the  compaction  gains  of  both 
the  linear  programming  and  the  window  method  versus  the  filter  order.  The  number  of  frequencies  used  m 
the  linear  programming  method  is  L  =  512  while  the  periodicity  used  in  the  window  method  is  L  =  2 N.  The 
windows  used  in  the  linear  programming  are  triangular  windows  with  symmetric  order  L-N- 1.  In  the  window 
method,  the  autocorrelation  sequence  is  first  windowed  by  a  triangular  window  of  symmetric  order  N  to  find 
fi{n)  and  then  the  window  is  reoptimized. 

From  the  figure  we  observe  that  if  the  order  is  high,  one  has  slightly  better  compaction  gains  using  the  window 
method.  This  implies  that,  if  one  optimizes  the  window,  there  is  no  need  to  use  large  number  of  frequencies  in 
the  linear  programming  method!  More  importantly,  there  is  no  need  to  use  the  linear  programming  technique  for 
high  filter  orders.  However,  it  should  be  emphasized  that  if  the  windows  are  optimized  in  the  linear  programming 
method,  one  can  get  slightly  better  compaction  gains  than  the  window  method.  In  Fig.  10(b),  we  show  the  plots 
of  the  compaction  gains  of  the  two  methods  for  various  values  of  M  for  a  fixed  filter  order  of  65.  We  observe  that 
the  window  method  performs  very  close  to  the  linear  programming  method  especially  for  low  values  of  M.  We 
show  the  upper  bounds  on  compaction  gains  in  both  plots.  The  upper  bound  in  the  first  plot  is  achieved  by  an 
ideal  compaction  filter  and  that  in  the  second  plot  is  achieved  by  a  maximal  eigenfilter  as  discussed  in  Sec.  1.2. 

IV.  PROPERTIES  OF  FIR  COMPACTION  FILTERS  AND  BOUNDS  ON  COMPACTION  GAINS 

In  this  section  we  will  give  a  number  of  results  pertaining  to  the  properties  of  optimum  FIR  compaction 
filters  and  the  corresponding  compaction  gains  for  a  given  filter  order  N  and  the  number  of  channels  M. 

1.  A  necessary  condition  on  the  compaction  filter  for  optimality.  For  an  FIR  compaction  filter  H(z) 
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to  be  optimum  it  is  necessary  that  the  Fourier  transform  of  the  sequence  r(n)g*(n)  attains  a  maximum  at  the 
frequency  w  =  0,  where  g(n)  is  the  impulse  response  of  the  product  filter  of  H(z).  To  see  this,  consider  the 
Fourier  transform  of  the  product  r(n)g* (n)  at  a  frequency  uq: 


r  Sxx(eju)G(e^~^)^ 

J —IT  “  ' 


(39) 


But  the  LHS  is  the  output  variance  of  a  valid  compaction  filter  whose  product  filter  is  g(n)eJUJ°n.  Since  g{n)  is 
optimum  for  the  given  autocorrelation  sequence  r(n),  it  follows  that  the  above  RHS  attains  a  maximum  when 
w0  =  0.  If  the  process  is  real,  then  by  considering  g(n)  =  g(n)  cos  u0n  one  can  arrive  at  the  same  result  with 


real  coefficient  filters. 

2.  Class  of  random  processes  that  have  the  same  optimum  compaction  filters  for  a  pair  (M,N). 
Let  us  consider  the  objective  in  the  time  domain: 


cr^=r(0)+  r(n)9m(n)  (40) 

njtkM 

From  this,  one  can  deduce  that  an  optimal  compaction  filter  does  not  depend  on  r(fcM)  and  it  continues  to  be 
optimal  for  a  modified  autocorrelation  sequence  of  the  form  f  (n)  =  c  r(n),  n  ^  kM,  where  c  >  0.  We  will  see 
later  (Sec.  V)  that  in  the  two-channel  case,  if  c  <  0,  then  H(-z)  is  the  optimum  solution  for  f (n). 

3.  Monotonic  behaviour  of  the  optimum  FIR  compaction  gain.  Let  Gopt{M,N)  denote  the  optimum 
compaction  gain  for  a  given  number  of  channels  M  and  FIR  filter  order  N.  It  is  then  clear  that 

Gopt(kM,N)  <Gopt((k  +  l)M,N),  *  =  1,2,...,  and  Gopt(M,N)  <  Gopt(M,N  +  1).  (41) 


4.  Bounds  in  terms  of  eigenvalues.  Let  fL(n)  be  any  Nyquist(M)  with  nonnegative  Fourier  series  coeffi¬ 
cients.  Assume  L>  N.  Then, 

,  Amai|r(n)/2(n)jo  <  Gopt(M,  N)  <  \max{r(n)}^  .  (42) 

Here  the  notation  Amar{r(n)}^  stands  for  the  maximum  eigenvalue  of  the  Hermitian  Toeplitz  matrix  whose 
first  row  is  [r(0)  r(l)  . . .  r(N)].  The  inequality  on  the  left  follows  because  the  product  filter  g(n)  =  w(n) fL(n) 
achieves  that  bound  by  choosing  w(n)  as  the  product  filter  of  the  maximal  eigenfilter  of  the  Hermitian  Toeplitz 
matrix  formed  by  the  sequence  r(n )/£(n).  For  the  inequality  on  the  right,  note  that  there  exists  an  integer 
k>l  such  that  kM  >  N.  FYom  Sec.  1.2  we  have  Gopt(kM,  N )  =  Amal{r(n)  which  is  called  the  KLT  gain. 

From  (41),  it  follows  that  Gopt(M,N)  <  Gopt(kM , N)  =  Am0I{r(n)  . 

If  we  replace  fL(n)  by  a  positive  definite  Nyquist (M)  sequence  f(n)  of  order  N,  the  inequality  on  the  left 
continues  to  be  valid  because  w(n)f(n)  is  still  a  product  filter  of  a  valid  compaction  filter.  To  see  this,  note 
that  the  sequence  /(n)  can  be  extended  to  an  infinite  sequence  (e.g.,  using  autoregressive  extrapolation)  such 
that  its  Fourier  transform  is  nonnegative.  Hence  the  product  w(n)f(n)  has  nonnegative  Fourier  transform.  The 
Nyquist(M)  property  of  the  product  follows  from  that  of  /(n). 
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5.  Upper  bound  by  M.  For  all  FIR  compaction  filters  we  have 

Gopt(M,N)<M 


(43) 


with  strict  inequality  as  long  as  Sxx(ejw)  is  not  a  line-spectral  process.  To  see  this,  first  observe  that  G{e]u  )<M. 
Hence,  G(eP)Ssx{e**)£  <  M  flv  S«(e*-)&  =  The  equality  holds  if  and  only  if  G(e^)  =  M 

for  all  u  for  which  Sxx{e*a)  #  0.  If  Sxx(e^)  is  not  line-spectral,  this  requires  G(eju)  to  be  identically  zero 
for  some  region  of  frequency  which  is  impossible  since  the  order  is  assumed  to  be  finite.  For  M  —  2,  we 
will  derive  another  upper  bound  for  Gopt(2,N)  in  Sec.  5.1  (see  (53)).  This  bound  will  be  applicable  under 
certain  conditions  on  the  input  power  spectrum,  to  be  made  more  explicit  in  that  section.  Under  some  further 
conditions  that  we  will  describe,  the  stated  bound  can  in  fact  be  achieved. 

V.  ANALYTICAL  METHOD 

In  this  section  we  consider  the  special  case  of  two  channels  (M  =  2)  and  assume  that  the  input  x(n)  is  real 
so  that  the  compaction  filter  coefficients  h(n)  can  be  assumed  to  be  real.  For  this  two-channel  case  we  will  show 
that  the  optimal  product  filter  G(eju)  can  sometimes  be  obtained  using  an  an^tical  method  instead  of  going 
through  a  numerical  optimization  procedure.  We  will  also  present  a  number  u.  examples  which  demonstrate 
that  the  method  is  applicable  to  a  large  class  of  input  power  spectra.  Also  presented  are  examples  where  the 
analytical  method  can  be  shown  to  fail. 

The  analytical  method  is  motivated  by  the  fact  that,  under  some  conditions  to  be  explained,  the  objective 
function  (5),  which  can  be  written  as  crj)  =  f*r  G(ejlJ)Sxx(eju)^  can  be  conveniently  expressed  as  a  finite 
summation.  The  summation  involves  a  modified  polyhase  component  of  G(ejlJ)  at  a  discrete  set  of  frequencies 
w*  determined  by  the  psd  Sxx(eju).  This  will  allow  us  to  optimize  the  modified  polyhase  component,  and  hence 
G(z),  essentially  by  inspection. 

The  inspiration  for  our  work  in  this  section  comes  from  the  recent  contribution  by  Aas  et  al.  [21]  where 
the  Gaussian  quadrature  technique  is  cleverly  used  to  adress  the  problem  of  optimizing  the  frequency  responses 
of  filters.  Our  work  in  this  section  differs  in  a  number  of  respects.  First  we  do  not  use  Gaussian  quadrature, 
but  take  advantage  of  an  elegant  representation  for  positive  definite  sequences  which  results  from  the  theory 
of  line-spectral  processes.  Second,  we  take  into  account  the  knowledge  of  the  input  psd  in  the  optimization 
process.  Finally  we  present  several  design  examples  demonstrating  the  usefulness  as  well  as  the  limitations  of 
the  proposed  method. 

Suppose  the  product  filter  G{z)  =  E^L-at  SO1)*-"  is  expressed  in  the  traditional  polyphase  form  [13] 
_  E0(z2)  +  z_1  Ei(z2).  For  the  real  coefficient  case  we  have  g(n)  =  g(-n),  and  it  follows  that  the 
coefficients  of  the  FIR  filter  Ei(z)  have  the  symmetry  demonstrated  in  Fig.  11.  This  implies,  in  particular,  that 
Ei(z)  =  0  for  z  =  -1.  By  factoring  the  zero  at  z  =  -1  we  can  write  Ei(z)  =  (1  4-  z)Gi(z)  where  Gi(z)  is  FIR 
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(44) 


and  has  symmetric  real  coefficients.  We  therefore  obtain  the  modified  polyphase  representation 

G(z)  =  1  +  (z  +  z-1)Gi(z2),  i.e.,  G(e*w)  =  1  +  2cosWG1(e>'2“’) 

where  we  have  also  used  the  Nyquist(2)  property  of  G(z).  Since  Nyquist  condition  and  nonnegativity  of  G(e^) 
together  imply  0  <  G(eju)  <  2,  we  have  the  following  bound  on  the  modified  polyphase  component  G i(eJW): 

_ 1 - <Gi(eju>)< - 7  7TT,  -ir<v<  7T  (45) 

2cos(w/2)  “  V  ;  "  2cos(u;/2) 

Notice  that  G(z)  and  Gi(z)  can  be  determined  from  each  other  uniquely.  We  shall  express  the  output  variance 
ct2  in  terms  of  Gi(eju)  so  that  we  can  see  how  to  optimize  the  coefficients  of  Gi(z).  For  this,  write  the  input 
psd  in  the  traditional  polyphase  form  as  Sxx(z)  =  50(z2)  +  z^S^z2).  Then  <r2  can  be  simplified  into  the  form 
a\  =  r(0)  +  J^Gi(e^x(e^)^  where  ¥*(z)  =  (1  +  z_1)Si(z)  or  equivalently 

^(e^)  =  cos(w/2) (sxx(e^2))  -  Sxx(ej^~u/2)))  (46) 

Using  Parseval’s  relation  the  objective  can  be  written  as 

(N— 1)/2 

cr2=r(0)+  5Z  3i(n)V>*(n)  (47) 

n=-(JV-l)/2 

where  ipx(n)  is  the  inverse  transform  of  'J'x(z)  which  is  produced  below  explicitly  for  convenience. 

ipx(0)  =  2r(l),  i/>x(l)  =  r(l)  +  r(3),  . . . ,  V’i(^— )  =  r(N  ~  2)  +  r(Ar)-  (48) 

Note  that  ipx(n)  =  ipx{-n).  Our  aim  is  to  maximize  the  second  term  in  the  expression  (47)  for  fixed  ipx{n) 
(i.e.,  fixed  input)  by  choosing  gi(n)  under  the  constraint  (45)  and  the  usual  filter-order  constraint.  Under  the 
assumption  that  the  input-dependent  sequence  tpx(n)  is  positive  or  negative  definite  (see  Sec.  1.1  for  definition) 
we  will  show  how  this  can  be  done  analytically.  (The  significance  of  this  assumption  on  ^(n)  is  explained  in 
Sec.  5.3).  We  will  need  the  following  representation  theorem  for  positive  definite  sequences: 

Theorem:  Representation  of  positive  definite  sequences.  Given  a  positive  definite  sequence  of  m  +  1 
complex  numbers  n  =  0, . . . ,  rn },  there  exists  a  representation  of  the  form 

m 

4>{n)  =  Y^akejUkn,  n  =  0,...,m  (49) 

*=o 

where  a^’s  are  all  positive  and  aVs  are  all  distinct. 

Proof.  Let  P  be  the  (m  + 1)  x  (m  + 1)  Hermitian  Toeplitz  matrix  whose  first  row  is  4>T  =  [4>(0)  . .  •  4>(jn) ]. 

Consider  the  extension  of  P  into  a  singular  (m  +  2)  x  (m  +  2)  Hermitian  Toeplitz  matrix  P  such  that  its 
(m  +  1)  x  (m  +  1)  principal  submatrix  is  P.  This  extension  is  merely  augmenting  an  extra  element  <j>{m  4- 1) 
to  the  end  of  $  and  forming  the  corresponding  Hermitian  Toeplitz  matrix.  The  number  <f>(m  +  1)  is  chosen  to 
make  P  singular.  This  can  always  be  done  because  of  the  following  reason:  for  the  matrix  P,  one  can  run  the 
well-known  Levinson  recursion  procedure  [26]  to  obtain  the  optimal  mth  order  predictor  polynomial  Am(z).  If 
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one  now  considers  the  following  continuation  of  the  recursion  Pc(z )  =  Am{z)  +  cz  (  +  Mm(z)  with  |c|  —  1, 
then  this  corresponds  to  the  singular  predictor  polynomial  of  a  random  process  with  singular  autocorrelation 
matrix  P.  The  result  now  follows  from  a  well  established  fact  [26,  27]  that  states  that  a  WSS  process  is  line 
spectral  with  exactly  m  +  1  lines  if  and  only  if  its  (m  +  1)  x  (m  +  1)  autocorrelation  matrix  is  nonsingular  and 
(m  +  2)  x  (m  +  2)  autocorrelation  matrix  is  singular.  ■ 

Remarks.  It  is  clear  that  Pc(z)  is  also  the  minimal  eigenfilter  of  P .  The  zeros  of  Pc(z)  are  all  on  the  unit  circle 
and  distinct.  Let  {ejUk ,  k  =  0, . . . ,  m}  be  these  zeros.  The  distinct  frequencies  {u>k,  k  =  0, . . . ,  m}  are  referred 
to  as  the  line-spectral  frequencies.  The  representation  (49)  is  not  unique  because  of  the  nonuniqueness  of  the 


unit  magnitude  constant  c  in  the  proof. 

Real  case.  The  predictor  polynomial  Am(z)  and  the  constant  c  are  real.  Hence  we  have  two  cases:  c  =  ±1. 
The  case  c  =  1  leads  to  a  symmetric  singular  predictor  polynomial  Pi(z),  while  the  case  c  =  -1  leads  to 
an  antisymmetric  polynomial  P-i{z).  It  is  a  well-known  fact  that  the  distinct  unit-circle  zeros  of  these  two 
polynomials  are  interleaved  [22],  For  simplicity  assume  that  m  is  odd.  Then  P-i(z)  has  both  of  the  zeros  z  =  1 
and  z  =  -1  and  Pi (z)  has  none  of  them.  Using  Pi(z),  we  have  the  following  representation  for  a  real  positive 

definite  sequence  <j>(n ): 


(m— 1)/2 

<p(n)=  ^2  @k  cosw/tn,  n  =  0,  ...,m 


(50) 


k=0 


where  /?*’s  are  all  positive  and  u!k’s  are  all  distinct  and  different  from  0  and  7 r. 


5. 1.  Derivation  of  the  analytical  method 

Assume  for  simplicity  ( N  - 1)/2  is  odd  and  let  {^(n),  n  =  0, ....  (AT  -  l)/2}  be  positive  definite.  Applying 


the  real  form  of  the  representation  we  have 

(n-3)/4  n_1 

tpz(n)  =  ^2  Pk  coswfcn,  n  =  0,..., — - — 

k=0 

The  objective  (47)  can  therefore  be  written  as 

(N-3)/4  (JV-l)/2  (AT— 3)/4 

o-2  =  r(0)+  J2  ^  £  ffi(n)  coswfcn  =  r(0)  +  J2  ^Gi(e}u,k) 

k=0  n=— (JV— 1)/2 


(51) 


(52) 


From  (45),  the  output  variance  (52)  is  maximized  if  Gi(e]u,k)  =  2cos(L;72y>  fc  =  0, . . ., (AT  -  3)/4.  This  implies 
G(e?k>k! 2)  =  2,  and  by  Nyquist(2)  property  G(e^K~IJk^)  =0,  k  =  0, . . . ,  (N  —  3)/4.  Notice  that  these  zeros 
are  all  located  in  the  region  (tt/2 ,7r).  Since  0  <  G(e^)  <  2,  the  derivatives  of  G{e?“)  should  vanish  at  the 
above  frequencies.  Hence  we  should  have  GV'Ufc/2)  =  0,  k  =  0, . . . ,  (AT  -  3)/4.  In  view  of  (44),  this  in  turn 
implies  G[(ejUk)  =  4^^^,  k  =  0, . . . ,  {N  —  3)/4.  The  total  number  of  constraints  on  Gi(eju)  and  Gi(eJW) 
is  N±l,  Since  Gi (z)  =  EiHN-i)/2«i(n)r"  with  9i(n)  =  9i(~n),  it  is  determined  uniquely.  If  these  ^ 
constraints  are  satisfied,  the  solution  G(e^)  so  found  is  necessarily  nonnegative  in  the  frequency  region  [tt/2,  tt] 
(Appendix  B).  If  it  is  nonnegative  in  the  region  [0,  tt/2)  as  well,  then  it  is  the  optimum  compaction  filter  with 
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the  corresponding  compaction  gain 


Gopt(2,N)  =  l  + 


^(N- 3)/4 

2-jk= 0 _ 2cos(tJfe/2) 


r(0) 


(53) 


If  however,  turns  out  to  be  negative  at  some  frequency  in  [0,7r/2),  then  it  is  not  a  valid  solution  and  the 

above  RHS  is  only  an  upper-bound  for  Gopt(2,N). 

Uniqueness  of  the  solution.  Assume  that  G(e*")  obtained  by  the  method  is  nonnegative.  Then  it  is 
the  unique  solution!  To  see  this,  assume  there  is  another  optimal  product  filter  K(z).  Assume  Ki{z)  is  its 
modified  polyphase  component.  Then,  there  exists  a  frequency  w*  among  the  line-spectral  frequencies  such 
that  K\  (e?Wk )  <  2cos(Lt/2)  •  Hence  the  summation  (52)  for  Ki(e?u)  is  necessarily  less  than  that  for  G\  (e-’~ ) . 
resulting  in  contradiction.  Notice  that  it  is  the  product  filter  G{z)  that  is  unique,  not  the  compaction  filter  H(z). 
However,  if  the  zeros  of  H(z)  are  constrained  to  be  on  and  inside  the  unit-circle,  then  H(z)  is  unique  as  well. 


5.2.  Construction  of  optimal  G(z) 

Consider  the  following  factorization  of  G{z): 

G(z)  =  G0(z)Gi(z)  (54) 

where  G0(z)  contains  the  unit-circle  zeros  determined  by  the  above  procedure.  Hence  we  have 

N- 3 

G0(z)  =  ]I  (z  +  2cos(wjt/2)  +  z"1)2  (55) 

Using  the  Nyquist(2)  property,  it  is  possible  to  determine  Gi(z)  and  hence  G(z).  Let  <?o(n)  and  ffi(n)  be 
the  impulse  responses  of  G0(z)  and  Gi(z)  respectively.  The  product  (54)  in  z-domain  is  equivalent  to  the 
convolution  in  time  domain.  Using  the  convolution  matrix  and  taking  into  account  the  symmetries  we  get 

g  =  Agj  (56) 

where  the  vectors  g,gi  have  the  components  =  g(2n)'§in  =  ffi(n),  n  =  0, . . . ,  (Ar-l)/2,  and  A  is  obtained 
from  the  impulse  response  g0(n).  From  the  Nyquist(2)  property,  it  is  clear  that  g  =  [1  0  . . .  0]T.  Hence  gi  (n)  is 
determined  by  inverting  the  above  equation.  The  invertibility  of  the  matrix  A  follows  from  the  fact  that  Gi  (z) 
is  uniquely  defined  and  hence  G(z)  is  uniquely  defined.  Therefore  a  unique  solution  to  gi  must  exist,  implying 
that  A  is  nonsingular. 

Efficient  determination  of  Gq(z):  We  will  show  that  we  can  obtain  Gq{z)  from  the  singular  predictor 
polynomial  Px(z)  without  having  to  find  its  roots.  For  this,  let  us  write  Px(z)  explicitly: 

N-3  ILzl 

Pi(z)  =  cz~ ^  f[(z- e>‘)(z  -  e~jUk)  =  cz~ ^  JJ  (z  -  2  cos  uk  +  ^_1)  (57) 

k=o  k=0 

Now,  consider  the  upsampled  polynomial  P\{z 2).  This  can  be  written  in  the  form  P\ (z  )  =  Pq(z)Po(—z),  where 
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(58) 


pQ(z)  is  a  polynomial  in  z~l  of  order  with  all  its  zeros  in  the  left  half  plane.  To  be  explicit: 

TV -3 

Po(z)  =  n  (z  +  2 cos(wfc/2)  +  z-1) 

k= 0 

Hence  we  can  write  G0(z)  =  z^P02(z).  Therefore,  given  the  singular  predictor  polynomial  Pi(z),  one  can 
apply  a  continuous-time  spectral  factorization  algorithm  [28]  to  Pi(z2)  to  obtain  Pq{z)  and  therefore  G0(z). 
Since  G(z)  can  be  determined  from  Gq(z)  as  we  explained  before,  we  observe  that  there  is  no  need  to  find  the 

roots  of  Pi(z)\ 

Spectral  factorization.  To  find  the  compaction  filter  H(z),  we  need  to  spectrally  factorize  G(z).  It  is  clear 
that  we  can  write  H(z)  as 

H{x)  =  HotfHiiz)  (59) 

where  H0(z)  and  Hi(z)  are  the  spectral  factors  of  G0(z)  and  Gi(z)  respectively.  We  can  deduce  H0{z)  immedi¬ 
ately:  Ho(z)  —  P0(z).  Hence  all  we  need  to  do  is  to  determine  H\ (z)  which  is  of  order  2  •  ^is  can  be  done 

by  a  discrete-time  spectral  factorization  of  Gi(z)  [29]. 

The  case  where  Pf1  is  even  can  be  treated  in  a  very  similar  manner.  In  this  case,  we  use  the  singular 
polynomial  P-i(z)  corresponding  to  c  =  —1  and  one  of  the  line-spectral  frequencies  is  0,  that  is,  z  —  1  is  a  root 
of  P-i(z).  The  resulting  product  filter  G(ejuJ)  continues  to  be  nonnegative  in  [tt/2,  tt].  We  skip  the  details  and 
give  the  summary  of  the  algorithm  for  both  cases. 

Summary  of  the  analytical  method 

Given  the  autocorrelation  sequence  r(n),  n  =  0,...,N,  where  N  is  odd,  we  would  like  to  find  an  optimum 
compaction  filter  H{z)  of  order  N.  We  first  obtain  the  sequence  ipx(n),  n  =  0,  ....(TV  — 1)/2  using  the  relations 
(48).  If  this  sequence  is  positive  definite,  then  we  do  the  following 

1.  Calculate  A*-i  (z).  the  optimum  predictor  polynomial  of  order  corresponding  to  the  positive  definite 

2 

sequence  and  obtain  Pc(z)  from 

Pc(z)  =  An=±(z)  +  cz-^A^iz-1)  (60) 

where  c  =  1  if  is  odd,  and  c  =  -1  otherwise. 

2.  Obtain  the  spectral  factor,  P0(z)  of  Pc(z2)  using  a  continuous  time  spectral  factorization  algorithm  and 
determine  Go(z)  =  z  *  P02(z)- 

3.  Calculate  Gi(z)  using  (56)  and  find  its  spectral  factor  i?i(z). 

4.  The  optimum  compaction  filter  is  H(z)  =  Po(z)H i  (z). 

Case  where  rpx(n)  is  negative  definite.  From  our  developments  for  the  positive  definite  case,  and  using 
the  sequence  -■<!> x(n ),  it  can  be  proven  that  the  optimum  compaction  filter  is  H(z)  =  H{-z )  where  H(z)  is 
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the  optimum  compaction  filter  for  the  positive  definite  sequence  xj>x(n)  =  -ipx(n).  However,  it  is  easier  to 
see  this  directly  by  looking  at  the  objective  in  time  domain  (40).  First  note  that  ij>x(n)  corresponds  to  the 
autocorrelation  sequence  f(n)  =  -r(n),  n  ^  0.  Let  g(n)  and  g(n)  be  the  product  filter  impulse  responses 
for  H{z)  and  H(z)  respectively.  The  objective  is  then  to  maximize  Y,n=i  -3(nK(n)>  This  has  the  solution 
—g(n)  =  g(n),  n  #  0.  Hence  we  have  G{z)  =  G(-z)  and  therefore  H{z)  =  H(-z). 

Remark.  The  observation  made  above  clearly  extends  to  the  case  where  f(n)  =  c  r(n),  c  <  0,  as  scaling  does 
not  change  the  optimal  solution  as  stated  in  Sec.  IV. 

Example  6:  AR(1)  process.  Let  N  =  3,  and  r(n)  =  pn  where  0  <  p  <  1.  Then,  ipx{0)  =  2 p  and 
•0i(!)  =  p(l  +  P2)-  The  Hermitian  Toeplitz  matrix  corresponding  to  {ipx(n),  n  =  0, 1}  is 


which  is  positive  definite.  Hence  we  can  apply  the  above  algorithm.  Running  the  Levinson  recursion,  we 
have:  Ai(z)  =  1  -  ^z~l  and  using  c  =  1  we  have:  Pi(z)  =  1  -  (1  +  p2)z~l  +  z"2.  By  straightforward 
calculation  P0(z)  =  1  +  y/z  +  />2z_1  +z~2  from  which  it  follows  that  G0(z)  =  (z+y/z  +  p2  +z-1)2  and  Gi(z)  = 
i  (-  _  2\/3  +  p2  +  z~l).  For  all  values  of  p,  Gi(z)  does  not  have  unit-circle  zeros  and  therefore  it  is 

2(3+p2)3/2  v  v  r  J  — — 

nonnegative.  The  spectral  factor  of  Gi(z)  turns  out  to  be  (a+fez"1)  where  a=\]  \/3  +  p2  +  \/2  +  P2 

and  b  =  -  yV3  +  p1  -  ^2  +  p2.  Finally  we  obtain 

H(z)  =  Fo(z)^i(z)  =  ^(3^p2)3/4  (a  +  (6  +  a%/3  +  P2)2_1  +  (a  +  by/ 3  +  p2)z~2  +  6z"3)  (62) 

The  product  filter  is  G(z)  =  w^37T(-z3+3(2+p2)z+2(3+p2)3/2+3(2+p2)z-l-z-3),  with  the  corresponding 
optimum  compaction  gain  Gopt( 2,3)  =  1  +  For  comparison  purposes  we  have  also  designed  compaction 

filters  using  the  linear  programming  and  the  window  methods  for  several  values  of  p.  Table  1  shows  the  filter 
coefficients  and  compaction  gains  using  the  linear  programming  and  the  window  methods  together  with  those 
obtained  by  the  above  analytical  method. 

Linear  phase  constraint  is  a  loss  of  generality.  Notice  that  the  compaction  filters  obtained  by  the  analytical 
method  in  Example  6  cannot  have  linear  phase.  This  is  observed  by  considering  the  zeros  of  G(z)  which  consist 
of  a  single  reciprocal  pair  and  two  double  unit-circle  zeros.  For  H(z)  to  have  a  linear  phase,  the  multiplicity 
of  the  zeros  of  the  reciprocal  pair  should  be  two.  Since  the  solution  G(z)  is  unique,  we  conclude  that  the 
linear-phase  constraint  on  the  compaction  filter  H  (z)  is  a  loss  of  generality. 

Example  7:  MA(1)  process.  Let  N  =■  3,  r(0)  =  1,  r(l)  =  p  >  0,  and  r(n)  =0,  n  >  1.  The  sequence  rpz(n) 
is  therefore  V’i(O)  =  2 p,  tpx(l)  =  p.  The  corresponding  Hermitian  Toeplitz  matrix  is 


(63) 


which  is  positive  definite.  Hence  applying  the  algorithm,  Aj  (z)  —  1  —  \ z  1  and  therefore  using  c  1,  the 
singular  predictor  polynomial  is  Pj(z)  =  1  -  z_1  +  z“2.  The  continuous  spectral  factor  P0(z)  of  Pi(z2)  is 
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P0(z)  =  1  +  2cos(7T/6)z-1  +  z~2  =  1  +  n/3z->  +  z‘2  and  therefore  G0(z)  =  (z  +  V*  +  z  1)2.  Hence  G1(z)  is 
calculated  to  be  Gx(z)  =  -*§ (z  -2^3  +  z~l).  The  spectral  factor  Hx{z)  of  Gi(z)  is  found  to  be 

Hl{z)  =  3-3/42-i /*(yfvf+Vi  -  \]%{3-\f2z-1)  (64) 


Hence  the  compaction  filter  H(z)  is 

3-3/42-i/2  ( ^  +  y/2  +  (^3  +  v/6  -  ^v/3  -v/2)z-x  +  (y/vi+Vt  -  \]  3  -  \/6)z~2  -  V^-V^z"3) 

V  (65) 


The  product  filter  is  G(z)  =  ~*fz3  +  f  z  +  1  +  and  the  corresponding  optimum  compaction 

gain  is  Gopt( 2,3)  =  1  + 

Example  8:  MA(1)  process,  general  order  N.  We  will  analytically  find  the  optimum  compaction  filters 
of  arbitrary  (odd)  order  N  for  MA(1)  processes.  Note  from  Sec.  IV  that,  if  we  find  H(z)  corresponding  to  a 
MA(1)  process  with  p  >  0  then  this  is  the  solution  for  all  MA(1)  processes  with  p  >  0  since  the  autocorrelation 
sequences  are  the  scaled  versions  of  each  other  in  the  sense  explained  in  that  section.  From  the  arguments  given 
for  the  case  where  V>*(n)  is  negative  definite,  H(-z)  is  the  solution  for  all  MA(1)  processes  with  p  <  0. 

Now,  following  the  steps  of  the  algorithm  we  have 

Pc(z)  =  1  -  z"1  +  z-2  -  . .  -  +  (-l)^z-^  (66) 

If  is  odd,  then  the  zeros  of  P1(z)  are  e±i-,  wk  =  (2k -  1)^,  k  =  l,...,(A  +  l)/4.  Therefore  the  roots 
of  Po(z),  hence  the  unit-circle  zeros  of  the  optimum  compaction  filter  H(z)  are 


e±jUk,  fit  =  7T  -  (2fc  -  1)- 


7T 


*  =  1,  —  , 


N  +  1 


(67) 


'N  +  3’  4 

Similarly,  if  js  even,  the  roots  of  P-i(z)  are  0,  e±JlJfc,  Wk  =  2 k N+3 ,  k  =  —  l)/4.  Therefore  the 

roots  of  P0(z),  hence  the  unit-circle  zeros  of  the  optimum  compaction  filter  H(z)  are 


rr,  e±jn‘ ,  fifc  =  7r  —  2A; 


7T 


AT +  3’ 


A-  1 


(68) 


The  rest  of  the  procedure  involves  spectral  factorization  and  it  is  not  easy  to  see  what  ffi(z)  will  be  in  closed 
form.  However  we  note  that  the  algorithm  successfully  finds  the  optimum  compaction  filter  for  any  odd  order 
N.  Table  2  shows  the  compaction  filters  and  the  corresponding  compaction  gains  for  various  filter  orders. 
Example  9:  KLT.  If  N  =  1,  then  we  have  Aq(z)  =  1,  and  P-i(z)  =  1  -  z-1.  and  therefore  P0(z)  = 
j  +  z-i;  Gq{z)  =  z +  2  + z-1,  and  Gi(z)  =  Hence  the  optimum  compaction  filter  of  first  order  is  H(z)  = 
J_(l  +  z-i)  if  r(l)  >  0  and  it  is  H{z)  =  -Wl-z"1)  if  r(l)  <  0.  Notice  that  these  correspond  to  the  two-channel 
KLT  subband  coder  which  is  known  to  be  fixed.  The  corresponding  compaction  gain  is  Gopt( 2, 1)  —  1  +  r(0)  • 
It  should  be  noted  that  the  above  filters  and  the  corresponding  compaction  gains  are  optimal  for  any  power 


spectrum  and  for  any  number  of  channels.  Hence  we  have: 

Gopt(M,l)  =  l  +  M,  M  >2 


(69) 
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If  r(m)  is  maximum  of  all  r(n)  where  n  is  not  a  multiple  of  M,  then  one  can  achieve  the  compaction  gain  of 
1  +  fcj^l  by  using  the  filter  ^=(1  +  z~m)  if  r(m)  >  0  and  the  filter  ^=(1  -  z~m)  if  r(m)  <  0. 

Case  where  tpx{n)  is  semidefinite.  Assume  that  x/>x(n)  is  positive  semidefinite.  Then  there  exists  an  integer 
P  <  (N  -  l)/2  such  that  {ipx(n),  n  =  0,1,...,P}  is  positive  definite  and  {ipx (n),  n  =  0,1,...,P  +  1}  is  only 
positive  semidefinite.  Then  we  can  replace  ( N  -  l)/2  in  the  above  arguments  with  P  and  write  the  objective 
(52)  in  terms  of  P  4-  1  corresponding  line-spectral  frequencies.  This  enables  us  to  determine  a  product  filter  of 
symmetric  order  2P  4-  1  <  IV.  If  this  resulting  filter  is  nonnegative,  then  we  have  found  the  unique  minimum 
symmetric  order  product  filter  that  is  optimum  among  the  filters  of  symmetric  order  less  than  or  equal  to  N\ 
The  case  where  ipx{n)  is  negative  semidefinite  is  similar,  the  details  are  omitted. 

Example  10:  Case  where  ^x(n)  is  positive  semidefinite.  Let  N  =  3,  r(0)  =  1  and  r(l)  =  r(3)  —  p>  0. 
Then,  ipx (0)  =  ipx(l)  =  2 p.  The  associated  Toeplitz  matrix  is 


2  P 


1 

1 


1 

1 


(70) 


which  is  positive  semidefinite  and  singular.  The  number  P  is  0  in  this  case  and  the  objective  (52)  is  l+2pGi(ej0). 
By  letting  G\  (ej0 )  =  § ,  the  product  filter  G  (z)  of  symmetric  order  1  can  easily  be  seen  to  be  f  z  + 1 4-  f  z  _  1 ,  and  it 
is  readily  verified  that  G(eju)  >  0.  In  fact  this  is  the  KLT  solution  with  the  compaction  filter  H(z)  =  ^(1+z-1). 
The  corresponding  optimum  compaction  gain  is  1  4-  p.  No  3rd  order  solution  can  achieve  better  gain  than  this. 


5.3.  Characterization  of  processes  for  which  the  analytical  method  is  applicable. 

For  the  technique  of  Sec.  5.1.  to  be  applicable  for  a  particular  order  N,  we  need  to  have  the  sequence  ipx(n) 
to  be  positive  or  negative  definite  for  that  N.  However,  al though  the  sequence  Vz(n)  is  positive  or  negative 
definite  and  the  algorithm  is  applicable,  it  may  happen  that  the  resulting  G(eJU)  is  not  nonnegative.  This  will 
be  the  case  if  G\  (eju),  obtained  from  nonnegative  Go(e’“),  is  not  nonnegative.  In  what  follows  we  will  describe 
a  class  of  processes  for  which  the  sequence  V;z  (n)  is  positive  definite  for  all  N  and  therefore  the  analytical 
method  is  applicable  for  all  orders  N. 

The  sequence  xpx(n)  is  positive  definite  for  all  N ,  if  and  only  if  'L x  (G^' )  is  not  a  line-spectrum  and  ^I(eJ  )  ^  0. 
Using  (46),  this  is  true  if  and  only  if  Sxx(eju)  is  not  a  line-spectrum  and 

(Sxx(ej“)  -  Sxx(ej^~u)))  >  0,  u  £  [0,  tt/2]  (71) 

We  will  say  that  the  process  is  ‘low-pass’  if  its  psd  satifies  the  above  condition.  Notice  that,  in  the  ideal  case, 
the  optimum  compaction  filter  for  that  type  of  process  is  the  ideal  half-band  low-pass  filter  [13]  (see  Fig.  12(a)). 
For  the  case  where  ipx(n)  is  negative  definite  for  all  N,  the  preceding  is  replaced  with 

(s**(eiM)  -  Sxxie^-^))  <0,u£  [0,  tt/2]  (72) 

For  this  type  of  input  process,  the  optimum  ideal  compaction  filter  is  the  ideal  half-band  high-pass  filter  [13]. 
Hence  we  will  say  that  the  process  is  ‘high-pass’  if  its  psd  satisfies  (72)  (see  Fig.  12(b)). 
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Cases  where  the  algorithm  fails.  Assume  that  the  process  is  such  that  the  sequence  {ipx(n),  n  0,...,(N 
1) /2}  is  positive  definite  and  therefore  the  algorithm  is  applicable  for  the  filter  order  N.  Assume  however  that 
one  of  the  line-spectral  frequencies  w*  is  close  to  n.  The  algorithm  will  require  e>(7r~a"‘/2)  to  be  a  zero  of  G(z). 
Hence  G{eju)  will  have  a  zero  close  to  tt/2.  But  this  may  be  impossible  if  the  order  N  is  low.  To  see  this, 
note  that  G(ejV/2)  =  1  from  the  Nyquist(2)  property  and  therefore  requiring  G(ejuJ)  to  have  a  zero  close  to  the 
frequency  7r/2  is  the  same  as  requiring  a  narrow  transition  band  for  G(eju)  which  is  impossible  if  the  order  is 
not  sufficiently  high.  One  can  however,  increase  the  filter  order  to  overcome  the  problem. 

Example  11.  Let  N  =  3,  and  r(n)  =  coswin,  wi  G  [0,7r/2).  Hence  rl>x( 0)  =  2coswi,  rpx(l)  =  cos3wi  +  cos u>i , 
and  i>z{n)  is  positive  definite.  Using  the  procedure  of  Sec.  5.2,  we  find  G0{z)  =  (z  +  2cosu;i  +  z  1)2  from 
which  it  follows  that  Gi(z)  =  -  ^oswx  +  z~l).  This  has  a  unit-circle  zero  if  wi  G  (zr/3,7r/2) 

and  therefore  Gi(eju)  is  not  nonnegative.  Hence  the  algorithm  fails  if  the  impulse  is  within  n/6  neighborhood 
of  7r /2.  We  have  designed  optimum  compaction  filters  for  the  above  autocorrelation  sequence  using  the  linear 
programming  method  for  various  values  of  u>i .  We  have  observed  that  the  optimum  compaction  filters  agree 
with  the  above  analytical  solution  if  wi  G  (0,  tt/3].  For  the  complementary  case  of  Wi  G  (-7r/3,7r/2)  where  the 
analytical  method  fails  for  N  =  3,  linear  programming  yields  the  solution  G{z)  =  -\zz  +  1  -  \z~z,  regardless 
of  the  exact  value  of  ui.  The  factors  G0{z)  and  Gi(z)  of  G(z)  are  G0(z)  =  (z  +  2cos(tt/3)  +  z'1)2,  and 
A  (7\  _  l _ (z  -  4cos('7t/3)  +  z-1).  This  is  the  same  as  the  previous  solution  except  that  in  the 

16  cos3 (tt/3)  '  \  /  J  ' 

previous  solution  is  replaced  with  a  constant  value  equal  to  7r/3. 

As  another  example,  let  us  fix  u*  =  2tt/5  >  tt/3,  and  find  the  optimal  FIR  compaction  filter  of  order  5.  The 
corresponding  product  filter  is  <?(z)  =  \z*  +  1  +  \z~z  and  the  compaction  gain  is  Gopt( 2, 5)  =  2  which  is  the 
largest  possible  gain  for  M  =  2!  Since  the  process  is  line-spectral,  this  is  not  surprising.  The  important  point 
here  is  that  while  the  algorithm  is  not  successful  for  the  filter  order  3,  it  is  successful  for  a  higher  order  5. 
Example  12:  Case  where  the  process  is  multiband.  Finally  we  will  consider  an  example  in  which  the 
input  is  not  “low-pass”  or  “high-pass”  but  rather  is  of  multiband  nature.  Let  r(0)  =  l,r(l)  =  ^,r(2)  =  0,  and 
r(3)  =  _1.  The  sequence  rpx(n)  is  positive  definite  for  N  =  3  so  that  the  algorithm  is  applicable.  There  is 
more  than  one  way  to  extrapolate  this  sequence  and  find  the  corresponding  psd.  For  example,  one  can  consider 
MA(3),  AR(3),  or  line-spectra(4).  In  all  three  cases,  we  have  verified  that  the  psd  is  neither  ‘low-pass  nor 
‘high-pass’.  Rather  it  is  of  multi-band  nature.  Applying  the  algorithm  steps  we  have  G0(z)  =  (z  +  75  +  z-1)2 
from  which  it  follows  that  <?i(z)  =  -\/2 (z  -  \/2  +  z"1).  This  has  a  complex  conjugate  pair  of  unit-circle 
zeros!  Hence  Gi(eju)  is  not  nonnegative  and  therefore  G(ejul)  =  Go(eiuJ)Gi(eju)  is  not  nonnegative  either.  The 
algorithm  halts  because  G\{e?tJ')  cannot  be  spectrally  factorized. 

VI.  CONCLUDING  REMARKS 

In  this  paper  we  have  developed  some  design  techniques  for  optimal  FIR  compaction  filters  and  elaborated 
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some  of  their  properties.  We  have  proposed  a  procedure  to  guarantee  the  nonnegativity  of  the  linear  program¬ 
ming  solutions.  Multistage  (IFIR)  designs  are  considered  and  some  design  examples  are  provided  to  demonstrate 
the  usefulness.  We  have  developed  a  new  and  efficient  design  method  that  we  called  the  window  method  which 
does  not  involve  any  of  the  iterative  optimization  techniques  such  as  linear  programming.  We  have  given  its 
relation  to  the  linear  programming  technique.  Finally  we  have  given  an  analytical  method  for  the  two-channel 
case  that  works  for  a  broad  class  of  random  processes.  In  all  the  techniques  we  have  concentrated  on  finding 
the  product  filter  G{z)  =  H(z)H(z)  of  the  compaction  filter  H(z).  Hence  the  final  stage  of  every  algorithm  is  a 
spectral  factorization  to  find  the  compaction  filter  H(z).  However,  we  can  first  calculate  the  compaction  gains 
for  different  filter  orders  using  the  product  filter  coefficients  g(n).  This  enables  one  to  decide  what  filter  order 
to  use.  The  bounds  we  gave  on  the  maximum  compaction  gain  are  also  useful  for  this  purpose. 


APPENDIX  A 


Fact  Al.  Let  xl{tl)  be  a  periodic  sequence  with  period  L  =  KM  and  consider  the  decimated  sequence 
yx(n)  =  xl(Mti)  which  is  periodic  with  period  K.  Let  Xi(k)  and  Fk(fc)  be  respective  Fourier  series 
coefficients.  Then  we  have  Yk(fc)  =  jj  ^2iLo 1  Xi(k  -F  iK). 

Proof.  The  Fourier  series  coefficients  of  yic(n)  is 


YK(k)  =  ^2  xL(Mn)W£n  =  5Z  7  S  X^(l)WL  MlnW^n  Xl ®  kYWk 

n=0  n=0  1=0  1=0  n=0 


(k-l)n 


1  M—l  K-l  1  ^v-x  1  * 

=  TjY.Il  Xrt + iK>w  7  wt‘” = 1 iK> 

■  i-0  7=0  n= 0  i=0 


n— 0  ^  l— 0 
K-l 


(73) 


»=0  j= 0 


M—l 

E 

i=0 


In  the  last  step  we  have  used  the  fact  that  J2n=o  Wj?"  —  " 

Proof  of  Lemma  2.  The  Fourier  series  coefficients  of  <5jr(n)  are  all  1.  Hence  from  Fact  Al,  it  follows  that 

2 =  xi{Mn)  —  if  and  only  if  Yx(k)  =  jj  ^i=o  xlQz  +  *-^0  ~  f-  * 


APPENDIX  B 

Proof  of  nonnegativity.  We  will  show  that  G(eJ'“)  obtained  by  the  procedure  in  Sec.  V  is  necessarily  non¬ 
negative  in  the  region  [ir/2,4  The  Nyquist(2)  property  of  G(e*w)  implies  G'(e^)  =  We  therefore 

have  G'(eJ’“*/2)  =  G'(eJ^-t'‘/2i)  =  0,  k  =  0,...,(N  -3)/4.  Now,  by  the  mean  value  theorem  in  calcu¬ 
lus,  we  also  have  G'(e^‘/2)  =  G'(ej^~^G))  =  0,  k  =  0,...,(N-  7)/4,  for  some  uk  e  (wfc,w*+i).  Notice 
that  since  s  are  aU  distinct  and  lie  in  the  open  region  (0,tt),  all  of  the  above  zeros  are  distinct.  The  total 
number  of  such  zeros  is  therefore  N  -  1.  Since  G(eju)  is  a  cosine  polynomial  of  order  N,  G'(e>w)  is  a  sine 
polynomial  of  order  N  and  therefore  it  can  be  written  in  the  form  G'(ejul)  =  sinu  T(cosw),  where  T(x)  is  a 
polynomial  of  order  N  -  1.  Excluding  the  zeros  at  0  and  7 r,  the  total  number  of  zeros  G'(e^)  can  have  in 
[0)7r]  is  iV  —  1.  Hence  G'(eju)  cannot  have  any  other  zero  on  the  unit-circle.  If  G{z)  has  a  zero  at  r  -  uk/2 
with  multiplicity  greater  than  2,  then,  G'(eju)  has  at  least  double  zero  at  that  frequency  implying  that  the 
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total  number  of  its  zeros  is  more  than  N  -  1  which  is  a  contradiction.  If  G(z)  has  a  single  zero  in  the  region 
(tt/2,  tt)  which  is  different  from  all  wk’s,  then,  by  applying  the  mean  value  theorem  once  more,  G'(eJU)  has 
to  have  another  zero  which  is  again  a  contradiction.  Hence  we  have  proved  that  G(e’“)  has  double  zeros  at 
tt  -  w*/2,  k  =  0, . . . ,  (N  -  3)  /4,  and  that  it  does  not  have  any  other  unit-circle  zeros  in  [tt/2,  tt]  .  This  in  partic¬ 
ular  implies  G(eju)  >  0  for  u>  €  [7r/2,7r].  The  proof  for  the  case  of  even  is  similar,  the  details  are  omitted.  ■ 
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Lisf  of  Figures 


Figure  1:  Pertaining  to  the  FIR  energy  compaction  problem. 

Figure  2:  Two-channel  orthonormal  filter  bank.  The  filter  H0(z)  determines  all  the  other  filters. 

Figure  3:  Windowing  of  the  linear  programming  solution.  The  window  w(n)  has  nonnegative  Fourier  transform. 
For  illustration  purposes  the  window  is  chosen  to  be  triangular. 

Figure  4:  The  psd  of  an  AR(5)  process,  and  the  magnitude  square  of  an  optimal  compaction  filter  for  N  —  65 
and  M  =  2,  designed  by  linear  programming.  The  parameter  L  is  1024  and  a  triangular  window  is  used. 
Figure  5:  IFIR  compaction  filter  design,  (a)  Basic  configuration,  (b)  Equivalent  system. 

Figure  6:  Special  IFIR  design  configuration  where  Hq(z)  is  a  valid  compaction  filter  for  (Mo,  N0)  and  Hi  (z) 
is  a  valid  compaction  filter  for  (a)  H0(z)  is  fixed,  Hi(z)  is  optimum  compaction  filter  for  x0(n),  (b) 

Hi(z)  is  fixed,  H0{z )  is  optimum  compaction  filter  for  xi (n). 

Figure  7:  Decomposition  of  g(n)  as  w{n) f l(ti)  where  W(eJW)  >  0  and  Fi{k )  >  0. 

Figure  8:  The  procedure  to  find  FL(k):  SL( 0)  is  maximum  among  {SL(iK)},  hence  FL( 0)  =  M,  FL(IK )  = 
0,1  ^  0.  5x,(l  +  .f0  is  maximum  among  (Sr,(l-f  iK)},  hence  Fl{1  +  K)  =  M,  Fl{  1  +  iK)  —  0,1^1,  and  so  on. 
Figure  9:  Periodical  expansion  of  r(n).  (a)  L  >  2N ,  (b)  L  =  2 N. 

Figure  10:  Comparison  of  the  window  and  linear  programming  methods.  The  input  power  spectrum  is  as 
shown  in  Fig.  4.  (a)  Compaction  gain  versus  filter  order,  M  —  2,  (b)  Compaction  gain  versus  number  of 

channels,  N  =  65. 

Figure  11:  Coefficients  of  the  polyphase  component  £u(z)  of  the  product  filter  G{z).  Because  of  the  symmetry 
9{n)  =  3(_n)>  we  have  -£i(-1)  =  0. 

Figure  12:  Illustration  of  low-pass  and  high-pass  power  spectra,  (a)  A  low-pass  power  spectral  density  and 
the  corresponding  optimum  ideal  compaction  filter  response,  (b)  A  high-pass  power  spectral  density  and  the 
corresponding  optimum  ideal  compaction  filter  response. 
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Figure  6:  Special  IFIR  design  configuration  where  H0(z)  is  a  valid  compaction  filter  for  (M0,N0)  and  Hx{z) 
is  a  valid  compaction  filter  for  (Mx,Ni).  (a)  H0{z)  is  fixed,  ifi(z)  is  optimum  compaction  filter  for  x0(n),  (b) 
Hi(z)  is  fixed,  H0(z)  is  optimum  compaction  filter  for  x\ (n). 


Figure  7:  Decomposition  of  g{n)  as  u)(n)/i(n)  where  W(e]u)  >  0  and  FL(k)  >  0. 
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Figure  8:  The  procedure  to  find  FL(k):  SL(0)  is  maximum  among  {SL(iF)},  hence  FL(0)  -  M,  FL{IK)  - 
0,1  ^0.  Sl(1  +  K)  is  maximum  among  {^(l  +  ii^T)},  hence  Fl(1  +  K)  =  M,  Fl{1  +  IK)  —0,1^1,  and  so  on. 
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Figure  9:  Periodical  expansion  of  r(n).  (a)  L  >  2N ,  (b)  L  —  2 N. 
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Figure  10:  Comparison  of  the  window  and  linear  programming  methods.  The  input  power  spectrum  is  as  shown 
in  Fig.  4.  (a)  Compaction  gain  versus  filter  order,  M  =  2,  (b)  Compaction  gain  versus  number  of  channels, 

N  =  65. 
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Figure  11:  Coefficients  of  the  polyphase  component  E\(z)  of  the  product  filter  G(z).  Because  of  the  symmetry 
g(n )  =  g{—n ),  we  have  Ei(-l)  =  0. 
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Figure  12:  Illustration  of  low-pass  and  high-pass  power  spectra,  (a)  A  low-pass  power  spectral  density  and 
the  corresponding  optimum  ideal  compaction  filter  response,  (b)  A  high-pass  power  spectral  density  and  the 
corresponding  optimum  ideal  compaction  filter  response. 


p= 

0.1 

n  1 

window  method 

linear  programming 

0.5494144350 

0.6940928372 

0.5839818982 

i 

0.7789293967 

0.7136056607 

0.7658293099 

2 

0.2470689810 

0.0680132766 

0.2140666953 

3 

-0.1742690225 

-0.0661535225 

-0.1632362113 

compaction  gain 

1.1153 

1.1078 

1.1151 

P  = 

0.5 

n 

analytical  method 

window  method 

linear  programming 

0.5308991349 

0.6817974052 

0.5693221037 

1 

0.7963487023 

0.7258587819 

0.7821354608 

2 

0.2411149862 

0.0663296736 

0.2047520302 

3 

-0.1607433241 

-0.0623033026 

-0.1490404954 

compaction  gain 

1.5547 

1.5283 

1.5537 

PS 

=  0.9 

1  n 

window  method 

linear  programming 

0.4938994371 

0.6550553981 

0.5605331011 

i 

0.8279263239 

0.7510864372 

0.8017336546 

2 

0.2281902949 

0.0620169861 

0.1699982390 

3 

-0.1361269173 

-0.0540877314 

-0.1188544843 

compaction  gain 

1.9222 

1.9118 

1.9207 

Table  1:  Compaction  filter  coefficients  h(n)  and  corresponding  gains  for  AR(1)  processes  with  various  p  values. 
The  filter  order  is  N  =  3  and  the  number  of  channels  is  M  =  2. 


n 


n 

N  =  3 

N  =  9 

N  =  15 

N  =  21 

0 

0.5502267080 

0.3472380509 

0.2619442448 

0.2135948251 

1 

0.7781380728 

0.7212669193 

0.6444985282 

0.5837095513 

2 

0.2473212614 

0.5313628729 

0.6197371546 

0.6522949264 

3 

-0.1748825411 

-0.0301418144 

0.1178983343 

0.2237822545 

4 

-0.2357012104 

-0.2498547909 

-0.2185003959 

5 

0.0008621669 

-0.1127984531 

-0.1849242462 

6 

0.1250275869 

0.1367316336 

0.1009864921 

7 

-0.0141611881 

0.0800586128 

0.1357184335 

8 

-0.0608205190 

-0.0879123348 

-0.0574793351 

9 

0.0292806975 

-0.0512638394 

-0.0996883819 

10 

0.0616834351 

0.0386787054 

11 

0.0272577935 

0.0732876341 

12 

-0.0441106561 

-0.0298852666 

13 

-0.0065141401 

-0.0529392251 

14 

0.0275486991 

0.0255464898 

15. 

-0.0111966480 

0.0362662187 

16 

-0.0230508869 

17 

-0.0216340483 

18 

0.0206081274 

19 

0.0077883397 

20 

-0.0156868986 

21 

0.0057402527 

p  =  0.1 

1.1155 

1.1244 

1.1260 

1.1266 

compaction  gains  p  =  0.3 

.  1.3464 

1.3732 

1.3798 

p  =  0.5 

1.5774 

1  1.6301 

1.6330 

Table  2:  Compaction  filter  coefficients  and  corresponding  gains  for  MA(1)  processes.  The  number  of  channels 
is  2. 
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